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1.  INTRODUCTION 
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This  report  is  a  incomplete  draft  of  a  paper  describing  the  physics, 
algorithms,  and  numerical  methods  used  in  the  two-dimensional  magnetohydro¬ 
dynamic  simulation  code  MACH2.  It  is  being  released  in  this  form  so  that 
any  aid  it  might  provide  to  users  of  the  code  will  be  available  as  soon  as 
possible.  The  author  hopes  that  the  final  version  of  this  report  will  be 
improved  by  this  early  review,  and  will  gladly  receive  any  comments  or 
reports  of  errors  which  the  reader  might  care  to  send  to  him. 

The  proposed  outline  for  missing  sections  of  the  report  is  as  follows: 


Section 


10.  ITERATION  PROCEDURES 

10.1  JACOBI  ITERATION 

10.2  MULTIGRIO  ACCELERATION 

10.3  THERMAL  DIFFUSION 

10.4  MAGNETIC  DIFFUSION 

10.5  CORRECTION  OF  NON-SOLENOIOAL  DEVIATION 

10.6  LAGRANGIAN  HYDRODYNAMICS 

11.  IMPLEMENTATION 

11.1  DATA  STRUCTURE 

11.1.1  Geometry  Description 

11.1.2  Spatially  Dependent  Physical  Quantities 

11.2  BOUNDARY  CONDITIONS 

11.2.1  Boundary  Condition  Control 

11.2.2  Application  by  Stride 

11.2.3  Use  of  BCPOINT  and  BCPNTRS 

11.3  DATA  FLOW  BY  PHYSICAL  PROCESS 


In  addition,  material  will  be  added  to  the  sections  dealing  with  time  and 
space  differencing  to  describe  the  recently  added  code  which  ensures  that 
the  poloidal  magnetic  field  remains  divergence  free.  In  addition. 
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appendices  connecting  the  theoretical  descriptions  of  this  report  with  the 
actual  subroutines  and  variaoles  of  the  code  will  be  included  in  the  final 
version  of  this  report. 


3 
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Even  in  this  draft,  the  author  wishes  to  acknowledge  the  dependence  of 
the  development  of  MACH2  on  the  contributions  and  strong  efforts  of  many 
others.  Without  the  technical  foundation  and  help  provided  by  Jerry 
Brackbill,  this  effort  could  not  even  have  been  attempted.  He  provided  the 
source  code  for  MOQUI,  out  of  which  MACH2  grew.  The  faith  and  support  of 
Bob  Reinovsky,  Jim  Degnan,  and  Bill  Baker  was  as  inspirational  as  it  was 
essential.  Moreover,  considering  how  seldom  experimentalists  are  willing 
to  support  applied  theoretical  efforts  such  as  this,  it  was  extraordinary. 
Their  foresight  and  wisdom  in  this  matter  is  a  clear  example  for  others  to 
emulate.  The  importance  of  the  technical  advice  and  moral  support  given  by 
Jim  Buff  and  Norm  Roderick  can  not  be  overstated.  Their  physical  insights 
and  experience  with  other  simulation  codes  often  provided  the  keys  to 
unravel  particularly  knotty  puzzles.  Peter  Turchi  played  a  special  role 
which  must  be  mentioned.  When  the  development  of  MACH2  entered  that  most 
difficult  stage  where  the  code  ran,  but  its  answers  were  not  yet  trusted, 
his  insight  and  understanding  of  magnetohydrodynamics  were  an  essential 
part  of  the  process  of  persuasion.  In  the  end,  the  success  of  the  develop¬ 
ment  effort  hinged  on  the  skills  and  efforts  of  those  who  actually  wrote 
and  debugged  code.  Bob  Peterkin  deserves  most  of  the  credit  for  the 
physics  code  while  to  Tony  Giancola  goes  credit  for  the  graphics,  restart, 
and  system  interface  code.  Both  made  significant  and  substantial  contri¬ 
butions  in  all  phases  of  the  development  from  the  design  stage  through 
initial  applications.  The  author  is  proud  to  be  associated  with  this  most 
capable  team. 
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2.  PHYSICAL  MODEL 


The  physical  models  in  the  code  today  are  largely  those  that  are 
required  to  generate  believable  solutions  to  fast  imploding  liner  or  plasma 
flow  switch  problems.  Specifically,  the  code  includes  ideal  MHD,  resistive 
diffusion  of  magnetic  field  and  consequent  joule  heating,  thermal  dif¬ 
fusion,  and  radiation  cooling.  The  differential  equations  expressing  the 
conservation  laws  relevant  to  these  processes  will  be  described  in  detail 
below.  In  what  follows  d/dt  represents  the  convective  derivative. 

2.1  THE  MHO  ELECTROMAGNETIC  FIELD  EQUATIONS 

As  explained  in  Jackson  (Ref.  1),  plasma  phenomena  with  characteristic 
times  much  longer  than  the  plasma  oscillation  period  can  be  analyzed  by 
neglecting  the  displacement  current;  this  is  the  MHD  approximation.  In 
rationalized  MKSA  units  then,  the  MHD  fields  are  described  by 

^  X  F  +  —  =0 

at 

V  X  B  =  Ug  J  (1) 

supplemented  by  an  equation  relating  J  and  E  that  depends  on  the  properties 
of  the  medium,  i.e..  Ohm's  law 


E  =  nJ  ”  V  X  B 

The  electromagnetic  units  used  in  MACH2  are  selected  so  that  Ug  =  1  ^nd  are 
related  to  the  MKSA  units  by  the  following  equations: 
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■"^“V^^oMKSA  ""^SA 


'^M2'V%MKSA  ‘^MKSA 


’M2 


’M2 


'oMKSA 


’MKSA 


V  '^oMKSA 


(2) 


In  this  system  of  units,  the  magnetic  energy  density  is 


D 

'  m 


(3) 


and  has  units  of  Joul es/meter^ .  Henceforth,  the  subscript  M2  will  be 
suppressed,  and  all  equations  written  in  these  units.  Thus,  the  field 
equations  in  MACH2  are 


«  E  .  II  =  0 
V  X  B"  =  J 

Application  of  Ohm's  law  for  the  moving  fluid  yields 


(4) 


If  =  -  7  X  [(n7  X  B)  -  (7x  B)]  (5) 

This  can  be  rewritten  as 

^  =  -  V  X  (nf  X  B)  -  (B  7  •  7  -  B  .  v7)  (6) 

where  the  left-hand  side  represents  the  field  transport,  the  first  term  of 
the  right-hand  side  represents  diffusion  of  field,  and  the  second,  the 
modification  of  field  intensity  by  the  divergence  of  the  velocity  field 
transverse  to 

The  resistivity  used  consists  of  two  terms;  one  is  classical  and  due 
to  particle-particle  interactions,  the  other  is  non-classical  (anomalous) 
and  intended  to  model  particle  interactions  with  the  turbulent  fields.  The 
classical  resistivity  may  be  taken  from  a  SESAME  format  table  or  computed 
from  a  Spitzer-like  formula;  the  anomalous  resistivity  is  always  computed 
from  a  formula  based  on  the  ion  acoustic  frequency. 


2.2  THE  MHD  MOMENTUM  EQUATION 

The  compressible  Navi er-Stokes  equations  supplemented  by  the  addition 
of  the  electromagnetic  force  are 


p  _  =  .vp  +  j  X  B  +  V  .  5 


It  is  possible,  and  convenient,  to  write  the  entire  right-hand  side  as  the 
divergence  of  a  tensor: 


(-Vp  +  J  X  B  -t-  V 


jif-lr-  [<-P  "»ij> 


-  (7  .  S), 


2.3  THE  CONTINUITY  EQUATION 


Conservation  of  mass  is  represented  by 


do  _  ~ 

-  -  -07  .  V 


2.4  THE  ENERGY  EQUATION 

If  we  let  e  be  the  internal  energy  per  unit  mass,  then  the  energy 
equation  in  MACH2  is 


0  ^  =  -p7  •  V  +  (nV  X  B)  •  V  X  B 


+  (a  •  V)  •  V  -  V  • 


and  includes  the  flow  work.  Joule  heating  due  to  diffusion  of  magnetic 
field,  diffusive  transport  of  energy  by  thermal  conduction  or  equilibrium 


diffusion  of  radiation,  and  radiation  cooling.  The  thermal  conduction  flux 

is  limited  to  some  fraction  of  the  free-streaming  flux  which  would  be 

achieved  if  all  the  plasma  on  the  high  temperature  side  moved  opposite  to 

the  gradient  at  the  local  average  thermal  velocity.  The  limit  on  the 

thermal  ^'lux  is  actually  applied  by  limiting  the  thermal  conductivity.  Let 

k  be  the  maximum  allowable  value  of  the  thermal  conductivity  determined 
max 


r  f 
Cl 

k  =  -  (ir 

max  |-^| 

Then,  if  the  thermal  conductivity  determined  by  local  plasma  conditions  is 
kp ,  the  effective,  or  limited,  thermal  conductivity  k^  is  determined  by 


"e  k  +  k 

p  max 

The  thermal  conductivity  k^  which  gives  an  energy  flux  equivalent  to  that 
due  to  diffusion  of  a  radiation  field  in  equilibrium  with  the  plasma  is 


The  diffusive  flux  of  energy  is  then 


F  k  +  k  1  vT 

diff  "  e  rx 


"he  radiation  cooling  rate  is  taken  to  be 


q  =  -C^  ^T 
^  rao  3 


This  is  applicable  when  the  plasma  is  thin,  and  the  rad-  tion  spectrum  may 
be  assumed  to  be  Planckian. 


2.5  THE  EQUATIONS  OF  STATE  AND  TRANSPORT  PROPERTIES 

All  of  the  equations  of  state  quantities,  and  some  of  the  transport 
properties,  are  determined  using  EOSPAC  (Ref.  2)  to  do  table  look-up  in  the 
Los  Alamos  SESAME  Equation  of  State  Library  (Ref.  3).  The  models  used  to 
generate  these  tables  are  among  the  best  known;  however,  there  are  regions 
in  density  and  temperature  space  where  interpolation  between  models  has 
been  used.  In  general,  those  regions  are  there  because  no  tractable 
numerical  model  exists  which  applies;  hence,  no  simple  analytic  model  is 
likely  to  be  superior! 


im.W/iSSS 
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3.  BOUNDARY  CONDITIONS 


Modeling  complex  experiments  requires  a  wide  variety  of  boundary 
treatments  even  for  this  modest  list  of  physics.  The  mathematical  expres¬ 
sion  of  these  will  be  detailed  below,  organized  by  the  physics  to  which 
they  apply.  Details  of  a  numerical  or  code  implementation  nature  will  be 
given  in  later  sections. 

The  philosophy  of  all  boundary  condition  application  in  the  code  is 
that  the  conditions  described  are  the  limit  of  conditions  in  the  fluid  as 
the  boundary  is  approached.  Thus,  the  boundary  of  the  computational  region 
is  the  edge  of  the  region  in  which  physical  modeling  is  done.  Quantities 
outside  this  region  are  assumed  known  or  related  to  quantities  inside  by, 
at  most,  simple  instantaneous  geometric  statements.  Specifically,  no  time 
integration  of  spatially  varying  quantities  is  carried  on  outside  the  prob¬ 
lem  boundaries. 

Figure  1  shows  a  small  section  of  problem  boundary  and  gives  the  geo¬ 
metric  quantities  which  will  be  used  in  the  description  of  boundary  con¬ 
ditions  to  follow. 


Figure  1.  Boundary  geometric  quantities. 
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The  tangent  vector  field,  t,  and  the  normal  vector  field,  n,  will  be  taken 


to  be  normalized  so  that 


Ini  =  F  =  1 


No  convention  will  be  adopted  for  the  sense  of  n  and  t;  it  will  be  speci' 
fied  only  if  circumstances  require  it. 


The  normal  and  tangent  vectors  always  are  taken  as  lying  in  the  compu 
tational  plane.  That  is 


n  .  e  =0 
9 


t  .  e  =  0 
9 


for  cylindrical  symmetry  and 


n  .  e  =0 
z 


t  •  e  =0 
z 


for  planar  symmetry. 


Discussion  in  the  remainder  of  this  section  will  be  restricted  to  the 


cylindrical  case.  In  all  cases,  the  results  of  that  discussion  may  be 
extended  to  planar  symmetry  by  the  transformation  of 


r  -*■  X 


z  >  y 


9  Z 


for  and  the  substitution  r  =  1  where  r  is  used  algebraically  and  not  as  a 
coordi nate. 
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3.1  MAGNETIC  FIELD  BOUNDARY  CONDITIONS 

The  boundary  conditions  used  in  the  code  represent  those  appropriate 
within  the  fluid  at  the  surface  of  idealized  conductors  and  insulators,  at 
a  symmetry  boundary,  or  at  the  axis  of  cylindrical  symmetry. 

3.1.1  Perfect  Conductor— As  shown  in  Jackson  (Ref.  4),  the  surface 
charge  and  current  instantly  nullify  any  non-zero  field  in  a  perfect  con¬ 
ductor.  Since  the  normal  component  of  B"  and  the  tangential  component  of  T 
must  be  continuous  through  the  surface,  it  follows  that  as  the  surface  is 
approached  from  outside  the  conductor  the  following  are  true  in  the  limit: 


n  •  B  =  0 


n  X 


r  =  0 


(20) 


Since  the  magnetic  field  is  the  fundamental  quantity  in  the  code  the  second 
equation  above  must  be  transformed  into  an  equation  involving  F.  The 
appropriate  equation  is 

n  X  E  =  n  X  (nJ  -  v  x  B)  =  0  (21) 

The  resistivity  will  be  considered  scalar  here.  Thus,  nx(nT)  =  n(^xF). 
Since  Fx(vxB)  *  (n  •  8)7- (n”  •  7)  B  and  TT  •  F  =  0  we  have 

n  (n  X  J)  +  (7  .  7)  F  =  0  (22) 

If  the  boundary  is  such  that  n  .  7  =  0,  i.e.,  it  is  either  a  free  slip  or  a 
no  slip  boundary,  then  this  becomes 


11 


Written  in  terms  of  B,  the  coordi nati zed  version  of  this  in  cylindrical 
coordinates  (where  3/39  =  0)  is 


!n  .  B  =  0 
and 

n”  •  V  (t  .  f)  =  0 
and 

Jr  n  .  f  (rBg)  =  0 


(24 


(25 


It  is  interesting  to  contemplate  what  effect  a  porous  perfect  conductor 
(one  for  which  n  •  v  ^  0)  would  have  on  the  field  adjacent  to  it,  but  in 
the  range  of  energy  densities  of  interest  to  us,  it  would  probably  plug  up 
immediately  with  ablated  mass.  Thus,  perfect  conductors  are  modeled  by 
Equations  23  through  25. 


3.1.2  Perfect  Insulator--No  current  flows  into  or  out  of  a  perfectly 
insulating  wall.  Thus,  "n  •  T  =  0  at  the  wall.  Since 


n  •  J  = 


where  t  •  e  =  0  and  t  is  tangent  to  the  wall,  it  is  clear  that  rB  is  con- 
0  0 
stant  along  an  insulating  wall.  The  constant,  of  course,  is  proportional 

to  the  total  current  flowing  through  the  (possibly  empty)  center  of  the 
(possibly)  annular  surface  of  revolution  forming  the  wall.  It  can,  there¬ 
fore,  vary  from  one  insulating  boundary  to  another. 


While  n  •  J  =  0  is  quite  restrictive  on  8,,  it  has  no  content  for 

9 

8^  and  8^  in  azimuthal ly  symmetric  geometry,  since  8^  and  B  affect  only 

J  and  ni  e  due  to  the  cylindrical  symmetry.  Thus,  the  appropriate  bound- 
0  -^0 

ary  condition  for  8^  and  8^  at  a  perfect  insulator  is  that  they  be  specif¬ 
ied  by  external  considerations.  For  example,  they  may  be  determined  from 
specified  current  flowing  in  conductors  such  as  Helmholtz  coils  outside  the 
region  of  interest.  Figure  2  shows  such  a  configuration. 


CONDUCTING 

FLUID 


SPECIFIED 


INSULATOR 

Figure  2.  A  problem  with  insulating  boundary  conditions. 

The  magnetic  field  at  boundary  points  such  as  P  may  be  determined  by 
applying  the  Biot-Savart  law  to  the  segmented  conducting  loops.  This  field 
will  diffuse  into  the  conducting  fluid,  and  induce  current  in  it. 


3.1.3  Symmetry  Surface  —  It  is  often  possible  to  take  advantage  of 
symmetry  in  the  problem  geometry  to  eliminate  half  the  problem.  Usually, 
this  is  done  by  applying  a  boundary  condition  which  implies  that  the  con¬ 
ditions  just  outside  the  problem  region  are  a  mirror  image  of  those  inside 
it.  For  the  magnetic  field,  the  appropriate  conditions  may  be  derived  from 
the  requirements  that  tangential  components  of  (J^,  J  )  and  (B^,  B^)  he 
zero  on  such  a  boundary,  and  the  normal  components  of  (B^,B  )  be  continu¬ 
ous.  The  first  of  these  requires  that 


i  n  .  V  (rBg)  =  0  (27) 

while  the  other  two  require 

t  .  B  =  0 

n  .  7  (n  .  B)  =  0  (28) 

where  we  make  the  assumption  then  n”  is  extended  off  the  boundary  so  then 
rT  •  ^  =  0. 

3.1.4  The  Axis  of  Cylindrical  Symmetry — At  r  =  0,  the  appropriate 
conditions  are  that 


B 
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=  0 


B  = 
r 


0 


J 


9 


=  0 


Si  nee 


38  3B 

j  =  _ 2 

9  3z  3r 


(29) 


(30) 


Equations  29  are  completely  equivalent  to 


3.2  HYDRODYNAMIC  BOUNDARY  CONDITIONS 


The  boundary  conditions  applied  to  the  hydrodynamic  quantities  repre¬ 
sent  those  conditions  within  the  fluid  at  the  surface  of  walls,  at  inlets, 
and  outlets  where  external  conditions  are  presumed  known,  at  inlets  and 
outlets  where  external  conditions  are  the  same  as  the  internal  conditions 
adjacent  to  the  surface,  and  at  the  axis  of  cylindrical  symmetry. 


3.2.1  Wal 1 s--At  a  boundary  representing  an  impermeable  wall,  the 
relevant  boundary  condition  is  either  free  slip 


n  •  V  =  0 


or  no  slip 


V  =  0 


3.2.2  Inlets  and  Outlets — At  an  inlet  the  properties  which  determine 
the  stress  tensor  and  the  convective  fluxes  must  be  specified.  There  the 
boundary  conditions  must  include 


e  =  e.  ,  .  or  T  =  T.  ,  ^ 
inlet  inlet 


inlet 


They  rust  also  include 


WvhI’ 


'tv 


The  same  quantities  must  be  rest^'icted  at  an  outlet; 


n  •  V  e  =  0 


n  •  7  0  ”  0 


n  .  7  V  =  0 


3.2.3  The  Axis  of  Cylindrical  Symmetry — At  r  =  0,  the  appropriate 


conditions  are 


V  =  0 

r 


'e  ■  ° 


3.3  THERMAL  CONDUCTION  BOUNDARY  CONDITIONS 

The  conditions  applied  to  the  temperature  at  boundaries  are  of  two 
types:  conduction  to  a  fixed  temperature  reservoir  of  infinite  heat 
capacity,  and  no  flux  at  all.  In  the  former  case,  the  condition  is 


T  =  T 


and  in  the  latter  it  is 


n  .  7  T  =  0 


4.  GEOMETRIC  PROBLEM  CLASS 


The  geometries  of  real  experimental  devices  are  usually  much  more  com¬ 
plex  than  the  idealized  geometries  that  motivate  the  experimental  design. 
Yet  these  deviations  from  the  ideal  geometry  may  have  significant  quantita¬ 
tive  and  qualitative  effects  on  the  results.  It  is  thus  essential  that  a 
two-dimensional  code  be  capable  of  solving  problems  in  domains  which 
closely  model  the  real  device.  Further,  some  interesting  experimental 
effects  cannot  be  seen  at  all  in  idealized  geometry.  Exploration  of  them 
with  simulation  codes  limited  to  domains  with  idealized  geometry  is  clearly 
impossi bl e. 

MACH2  was  designed  to  handle  any  of  a  broad  geometric  class  of  domains 
without  code  modification.  In  this  section,  that  geometric  class  will  be 
described,  and  some  complex  experimental  domains  shown  to  be  admissible  in 
that  class. 

A  subset  B  of  R'  given  by 

®a,b  "  I  ^  ^  %  5  y  5  (^0^ 


where  a  =  (a  , 

X 

to  as  a  bl ock . 
smooth  mapping. 


a  )  and  b  =  (b  ,  b  ),  is  a  rectangle  and  will  be  referred 

y  ^  y 

Recall  that  a  di ffeomorphism  is  a  smoothly  invertible 
that  is,  a  sufficiently  differentiable  function 


f  :  U  C 

for  whicn  there  exists  sufficiently  differentiable  function 

2  2 
g  :  W  C  R  >  R 

called  the  inverse  of  f,  which  satisfies 
g( f (w))  =  u  for  al 1  ueU 
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(41) 


(42) 


(43) 


f(  g(w))  =  w 


for  all  weW 


A  region  BCR  which  is  the  di f feomorphi c  image  of  a  block  will  be 
called  b1 ock-1 ike. 

The  regions  shown  in  Figure  3  are  examples  of  block-like  regions, 
which  may  be  thought  of  as  rectangles  formed  from  rubber  sheets  which  are 
deformed  by  the  di f feomorphi sm.  The  examples  in  Figure  A  are  also  block¬ 
like  except  that  the  di ffeomorphism  is  not  one-to-one  on  the  entire  block, 
but  is  multiple  valued  on  at  least  one  edge.  Thus,  the  inverse  mapping  is 
singular  at  the  sector  center  and  at  the  joining  cut  in  the  annulus.  Such 
regions  will  be  called  singular  block-like  regions. 

The  geometric  objects  which  MACH2  is  designed  to  treat  as  domains  con¬ 
sist  of  a  particular  kind  of  union  of  block-like  regions.  They  are  best 
described  as  (almost)  di ffeomorphic  images  of  a  union  of  blocks 


i  =  1 


with  the  a^  and  b. ,  chosen  so  that  any  two  blocks  which  meet,  do  so  only 
along  an  entire  edge  or  at  a  vertex.  Such  a  collection  of  blocks  will  be 
called  a  block  complex.  Figure  5  shows  some  examples  of  block  complexes. 
They  resemble  subsets  of  checkerboard  consisting  only  of  whole  squares.  It 
is  natural  to  call  the  di ffeomorphic  image  of  a  block  complex,  a  block-1  ike 
complex . 

This  class  is  quite  broad.  It  is  more  than  broad  enough  to  include 
realistic  experimental  configurations.  Figure  6  shows  two  configurations 
in  interest  which  are  within  the  geometric  class  described  above.  In  the 
present  version  MACH2,  corners  of  a  logical  block  must  meet  other  logical 
blocks  only  at  their  corners.  Further,  any  change  in  boundary  condition. 


DIFFEOMORPHISM;  X  (77,  ^  ) ,  Y 
(NUMERICALLY  GENERATED  GRID:  (X 


whether  physical  or  grid,  any  sharp  change  of  direction  of  the  region 
boundary,  or  any  change  in  uniform  initial  condition,  must  occur  at  a  block 
corner  or  edge.  Thus,  some  of  the  block  interfaces  shown  in  Figure  5  are 
required  for  non-geometric  reasons. 
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5.  COORDINATE  SYSTE^i 


In  MACH2  the  physical  quantities  are  computed  on  a  coordinate  system 
which  moves  arbitrarily  in  space.  Since  such  coordinate  systems  can  be 
made  to  remain  fixed  in  space  or  to  move  with  the  fluid,  they  are  called 
Arbitrary  Lagrangian-Eulerian  systems.  Clearly,  this  name  is  unfortunate, 
and  we  will  prefer  "arbitrary  coordinate  system"  or  "arbitrary  fluid 
descri pti on" . 

The  principal  advantages  of  such  a  description,  described  in  Refer¬ 
ences  5  and  6,  are  that  some  of  the  catastrophic  coordinate  distortion 
common  with  a  Lagrangian  description  and  the  diffusive  inaccuracy  of  the 
Eulerian  description  can  be  avoided.  The  principal  disadvantage  is  that  a 
prescription  for  the  coordinate  system  motion  must  be  supplied. 

In  MACH2,  the  coordinate  system  for  a  block-like  complex  is  described 
by  a  collection  of  function  pairs 

((x^(n,C,t)  ,  y^(n,C,t)}  |  (n.c)  e  B^,  t  e  R+,  1  =  1 . L}  (46) 

defined  on  the  collection  of  L  individual  blocks,  ,  making  up  the  block 
complex  which  is  the  domain.  The  coordinates  (n.c)  are  just  the  contin¬ 
uous  extensions  of  the  integer  array  coordinates  (i,j)  that  will  make  their 
appearance  upon  discretization.  All  of  the  physical  quantities  are  also 
described  as  functions  of  (n.c)  e  B^,  so  that  the  full  problem  is  actually 
given  by  a  parametric  description.  That  is  F,  p,  and  V  are  stored  as 
functions  of  1,  and  of  (n  ,C )  rather  than  (x,y). 

In  the  next  two  subsections,  the  range  of  choices  of  coordinate 
systems  available  in  MACH2  and  the  effects  of  this  choice  on  the  physical 
model,  as  represented  by  the  partial  differential  equations  of  Section  2, 
will  be  described.  Coordinates  in  the  logical  domain  (l,niC)  will  be 
represented  by  ^  . 


5.1  COORDINATE  SYSTEM  MOTION 


The  coordinate  system  x(5,t)  is  allowed  to  move  according  to 


^  (f  ,t)  =  a  (Xj(c,t)  -  x(c,t))  +  8v  (47) 

where  V  is  the  fluid  velocity,  is  an  "ideal"  coordinate  created  by  a 
numerical  grid  generation  procedure  which  we  will  describe  shortly,  and  the 
control  functions  aiX)  and  8 (F)  are  user  specified. 

The  Langrangian  coordinate  system  results  if  a  =  0  and  8  ?  1,  while 
the  Eulerian  coordinate  system  appears  if  both  a  and  8  are  identically 
zero.  Sensible  values  for  a  lie  in  the  range  0  <  a  <  1/T,  where  T  is  some 
relaxation  time,  while  8  must  be  either  0  or  1.  Thus,  a  is  a  relaxation 
parameter  telling  how  quickly  the  actual  coordinate  system  moves  toward  the 
ideal  coordinate  system. 

5.2  THE  IDEAL  COORDINATE  SYSTEM 


The  ideal  coordinate  system  is  generated  by  solving  a  Brackbill- 
Saltzman  variational  problem  on  the  block  complex  domain  subject  to  user 
specified  external  and  internal  boundary  conditions  and  using  a  weight 
function  selected  from  a  standard  family  by  the  values  of  user  specified 
parameters.  The  details  of  the  variational  problem  and  its  Euler  equations 
are  described  in  Brackbill  (Ref.  7),  and  their  implementation  for  block¬ 
like  complexes  is  described  in  Reference  8.  For  completeness,  a  brief 
description  will  be  given  here. 

Let 


/ax 

1 

a  X 

'  ac 

an 

1  IL 

iX, 

\ac 

an  1 
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so  that  J,  the  Jacobian  of  the  mapping  (n  .c  )  — ►  (x,y ) ,  is  given  by 


J  =  det  [ d^x) 

Then,  the  variational  problem  for  the  mesh  is  stated  as 


(49) 


Xo  ^  (d-Xj) 


dA  =  0 


(50' 


where  and  x^  are  user  specified  parameters,  D  is  the  block  complex 
domain,  and  the  nonlinear  operators  /  ,  and  are  local  measures  of  the 
mesh  smoothness,  its  nearness  to  desired  concentration,  and  its  orthogonal 
ity,  respectively.  These  operators  are  given  by 


//  (d-x)  = 


(— 1 


'^35'* 


f— 1 


(^) 

'^3n'' 


/J 


(d^x,x)  =  w(x,y) 


f  (d^xj 


fiA  iiL  +  iX  iXx 
'■35  3n  35  3n^ 


(51 


where  w  is  a  user  specified  volume  control  function.  Thus,  for  =  X^  =  ( 
the  solution  to  this  variational  problem  is,  in  some  sense,  the  smoothest 
coordinate  system  filling  the  physical  region;  whereas,  for  1  and 
x^  =  0  it  is  some  blend  of  a  smooth  coordinate  system  and  an  orthogonal 
one.  Of  course,  for  X^  0  the  character  of  w  will  have  an  effect  on  the 
solution  to  the  variational  problem,  so  if  w  depends  on  the  current  state 
of  the  simulation  variables,  then  the  coordinate  system  may  dynamically 
adapt  as  the  simulation  proceeds. 


The  role  of  the  function  w  in  controlling  the  adaptivity  of  the  mesh 
can  be  made  more  precise.  The  Euler  equations  for  Equation  49  with  x  = 


00 


+  2  [(wJ)  y  -  (wJ)  y  1  =  0 


I7  *  2  [x^(wJ)^-  •  0 


If  we  suppose  that  [x(y,5),  y(n.c))  satisfying  these  has  been  found,  let 

:/  (x,y)  =  X  (n.c)  y  (n  ,? )  -  x  (n,?)  y^(n,c).  and  transform  Equation  52 
C  n  n  C 

using 


_ 

3 

X  — 

+  V 

3 

C  3x 

3y 

X 

+  V 

3 

3n 

ri  3X 

3n 

then  we  obtain 

Thus,  the  Jacobian  of  the  coordinate  system  defined  by  Equation  52  must  be 
proportional  to  w  .  Since  the  Jacobian  of  the  coordinate  system  is 
approximately  the  cell  area  for  the  mesh,  cells  will  be  smaller  where  w  is 
large,  so  long  as  >  0. 

For  A^  of  order  unity,  the  terms  in  the  Euler  Equations  for  Equation 
50  which  are  derived  from  .7^  act  to  smooth  the  coordinate  system  and  reduce 
the  influence  of  the  weight  function.  The  resulting  coordinate  system  is 
thus  a  blend  of  those  determined  by  the  operators  in  Equation  51. 

5.3  THE  MATERIAL  DERIVATIVE 

The  only  terms  of  the  physical  model  of  Section  2  affected  by  the 
choice  of  coordinate  system  are  those  involving  the  material  derivative. 

The  material  derivative  for  a  quantity  in  Lagrangian  description  is 


^  ,  u. 

dt  3 1 


while  for  an  Eulerian  description  it  is 


For  the  arbitrary  fluid  description,  it  is 


3f, 

_  =  +  r  v-v  )  •  V  f 

dt  3 1  C''  a 


where  is  the  velocity  of  the  coordinate  system  defined  as  the  velocity 
of  points  in  physical  space  represented  with  fixed  coordinates.  Clearly, 
the  special  cases  ^  3nd  =  0  produce  Equation  55  and  Equation  56 
from  Equation  57.  Equation  57  is  derived  from  Equation  55  by  noting  that 


(n.C.t)  =  f  [x(n,c,t)  ,  y(n,C»t)  ,  t) 


and  differentiating  with  respect  to  t  to  obtain 


3f  3f  3f  3f 

d  ^  0  0  0 

at  ~  at  ^  \  a x  ^  ^t  ay 


Since  ( x^ ,  y^}  =  v^,  solving  Equation  59  for  af^/at  and  substituting  the 
result  into  Equation  55  produces  Equation  57. 

The  quantities  of  which  material  derivatives  are  required  in  the  model 
include  p,  e,  F,  and  7.  It  is  crucial  to  the  performance  of  the  code  that 
these  be  performed  in  a  way  which  conserves  certain  integral  quantities. 
The  details  of  the  spatial  differencing  by  which  that  is  accomplished  will 
be  covered  in  Section  7.  The  equivalent  integral  statements  of  the  differ 
ential  equation 


differenced  there  will  be  derived  here. 


If  R(t)  is  a  moving  region  of  space,  then 

+  7  •  ($ v)  =  0  (61 ) 

implies 

3T/R(t)  "  /3R(t)  ^(^s  -  V)  .  n  dA  (62) 

where  iji  is  described  in  Eulerian  coordinates,  7  is  the  fluid  velocity,  and 
v^  is  the  velocity  of  the  surface.  Conversely,  if  Equation  62  is  true  for 
all  regions  R(t)  then  Equation  61  is  a  consequence. 

To  illustrate  the  equivalence  of  these  two  statements,  we  must  param¬ 
eterize  the  moving  region  R(t)  by  x{T,t) 

7  :  Q  — ►R(t)  (63) 

where  x,(s,0)  =  s,  so  that  n  =  R(0).  Then 


^R(t)‘^^'^  ^  /q  p(y.t)  det  (d^)  ds 
in  which  d-y  is  the  Jacobian  matrix 


(64) 


^  s^  ^  i  j 


3yi 


so  that  det  (d-y)  is  the  Jacobian.  From  Equation  64 


(65) 


It  /R(t)  "  /q  n  ^  n  •  ^  ♦  It 
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—  [det  dry]  -  (7  .  v  )  (det  d^) 


lr/R(t)  "  Uit)  n  ^ 


Substituting  v  +  (v  -  v)  for  v  and  applying  3ij)/3t  +  7(4iv)  =  0  yields 


3T/R(t)  =  ^R(t)  ^  •  ^("s- 


An  application  of  Stokes  Theorem  then  produces  Equation  62.  Proof  of  the 
converse  uses  the  same  algebra  and  the  fact  that 

/j^  f  dV  =  0  for  all  R  (7C 

implies  that  f  nust  be  zero  if  it  is  continuous. 


By  applying  Equation  61  repeatedly  with 
be  shown  that  the  system 


=  p,  pe,  pv,  and  pB,  it  can 


p  V  •  V 


^  n 

dt  dt  " 


is  equivalent  to 


3t  '  R(t) 


0  eu  V 


■  3R(t) 


pc  V  V 


S 


V  } 


ri 


n  -'Rd)  '  o«(t)  O'  <*s 

It /R(t)  ■  !mt)  ''s  "  x*  ('2' 

It  is  clear  that  the  first  of  Equations  72  is  equivalent  to  the  first  of 
Equations  71.  The  second  equation  follows  from  expanding 

(pe)  +  V  •  (pev)  =  P  ^  ®  p7  •  v)  (73) 

The  second  of  Equations  72  implies  that  the  left-hand  side  of  this  is  zero; 
the  factor  in  parenthesis  on  the  right-hand  side  is  zero  by  the  previous 
remark.  Since  p  cannot  be  zero  anywhere,  the  second  of  Equations  71 
follows.  The  remainder  follow  in  a  similar  manner. 
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6.  TIME  DIFFERENCING 


The  physical  model  described  in  the  previous  sections  is  solved 
numerically  in  arbitrary  coordinates  by  a  time-split,  time-marching  algor¬ 
ithm,  The  physical  processes  which  are  time-split  are: 

Radiation  Cooling 

Thermal  and  Equilibrium  Radiation  Diffusion 

Resistive  Diffusion 

Lagrangian  Hydrodynamics 

Coordinate  System  Motion  and  Convective  Transport 

The  resistive  diffusion  and  the  Lagrangian  hydrodynamics  are  done  with 
implicit  time  differencing,  but  the  remainder  are  done  with  explicit  dif¬ 
ferencing.  The  time-step  is  selected  to  maintain  the  stability  of  the 
calculation. 

6.1  TIME-SPLITTING 

Often  known  as  fractional  time-stepping,  time-splitting  is  an  approach 
to  time  advancing  evolution  equations.  It  consists  of  the  sequential 
application  of  separate  portions  of  the  system  of  equations,  rather  than 
simultaneous  application  of  the  entire  set.  Thus,  the  conditions  or  state 
which  result  from  the  application  of  the  first  subset  of  equations  or  terms 
are  used  in  the  second,  to  compute  the  changes  in  state  dictated  by  that 
piece  or  pieces  of  the  governing  physical  model.  In  MACH2,  for  example, 
the  magnetic  field  which  results  from  the  action  of  the  resistive  diffusion 
for  the  current  time-step  is  used  for  the  initial  state  for  the  action  of 
the  Lagrangian  hydrodynamics  for  the  same  time-step.  The  inaccuracy  of 
this  time  differencing  can,  of  course,  be  made  small  by  making  the  time- 
step  sufficiently  small.  Later  subsections  of  this  section  will  describe 
precisely  how  the  various  equations,  and  in  fact,  even  terms  of  equations 
are  time-split. 


6.2  EXPLICIT  VERSUS  IMPLICIT  TIME  DIFFERENCING 


When  an  evolution  equation  of  the  form 


= 

9t 


is  solved  numerically,  the  values  of  the  dependent  variables  will  only  be 

computed  at  discrete  time  t^ ,  t2  , . .  .t|^ , . . . .  For  a  finite  difference  code, 

the  dependent  variable  x  is  a  vector  with  many  components,  one  for  each 

physical  variable  of  interest  (B  ,  8  ,  B  ,  p,  u,  v,  w  ,  x,  y,  e,  and  others 

r  9  2 

as  well  in  the  present  case)  at  each  spatial  location  of  interest  i.e.,  at 
each  grid  point.  If  we  use  to  stand  for  x(t  )  one  way  to 
difference  Equation  74  is 

i«„.l  -  ■  Ffx„)  (75) 

This  is  known  as  explicit  time  differencing,  because  the  new  time  value 
x^^^  is  given  explicitly  as  a  function  of  the  old  time  value  Xp  by 

X,=x‘*'*^tF7x)  (76) 

The  iterative  procedure  described  by  the  spatially  differenced  form  of  this 
equation  can  be  unstable.  That  is,  it  may  describe  non-physical  divergence 
from  physically  reasonable  behavior.  This  certainly  will  happen  if  the 
linear  part  of  Equation  76  amplifies  disturbances  with  characteristic 
lengths  of  variation  as  small  as  the  spatial  descreti zati on  dx.  This 
behavior  is  associated  with  the  existence  of  eigenvalues  of  the  lineari¬ 
zation  of  Equation  76  which  have  magnitude  greater  than  1.  Often  the 
requirement  that  these  eigenvalues  be  less  than  1  in  magnitude  can  be  main¬ 
tained  by  keeping  dt  small  enough.  The  existence  of  such  a  stability  limit 
on  the  time  step  is  typical  of  explicit  time  differencing;  implicit  time 
differencing  can  be  free  of  such  limits. 


Fully  it^plicit  time  differencing  for  Equation  54  is  given  by 


77' 


and  is  so-called  because  determined  implicitly  from  by  the  rela- 

t  i  on 


■  n+1 


-  dt  Ff, 


n+1^ 


(78 


This  must  be  solved  for  to  determine  the  new  time  values.  In  any 

interesting  case,  each  equation  represented  by  Equation  78  depends  on  many 
of  the  components  of  ot  which  there  may  be  hundreds.  Numerical 

solution  of  such  systems  is  typically  very  similar  to  the  solution  of  the 
steady-state  problem  (derived  from  Equation  77  by  taking  the  limit  dt  +  «) 

P[x^]  =  0  (79) 


If  the  partial  differential  equations  represented  by  Equation  74  are  para¬ 
bolic,  such  as  the  equation  of  thermal  diffusion,  then  the  steady-state 
problem  will  be  elliptic,  as  will  the  time  advance  problem.  Equation  77. 
Techniques  for  solving  elliptic  problems  numerically  such  as  successive 
over-relaxation,  alternating  direction  implicit  methods,  multigrid  methods, 
and  preconditioned  conjugate  gradient  methods  can  all  be  used.  Since  the 
operator  F  is  often  nonlinear,  these  methods  all  must  be  modified  to 
include  same  iterative  methods  for  the  solution  of  nonlinear  equations, 
usually  similar  to  Newton's  method. 

The  relationship  of  implicit  or  explicit  time  differencing  to  permis¬ 
sible  time-step  is  determined  by  stability  considerations  in  a  very  general 
framework.  In  either  case,  we  must  consider  the  linear  part  of  the  iter¬ 
ation  operator  L  such  that 


For  explicit  iteration  we  have  L=L  =l+dtdF  from  Equation  76;  the 

implicit  iteration,  Equation  73,  implies  L  =  L  h  7 1  -dt  d  F)”^.  Thus,  if 

“  '  X 

we  write  for  the  spectrum  o^  an  operator  M,  we  have  the  spectral  map- 

pi  ngs 


A  f  L^l  =  1  +  dt  Af  d  F) 
and 

'^i^  "  1  -  dt  A(d^F] 

For  parabolic  Equations  74,  A(d  F)  c  LHP  =  (ze  Cj  Rg(z)  <  01;  further 

if  d  F  is  taken  as  referring  to  the  discretized  F,  the  Afd  FI  is  bounded 
X  'X 

since  it  has  only  finitely  many  members.  Then  Equation  66a  implies  that  if 

dt  is  small  enough,  a(L^)  <  S(0,1).  Here  S[z^,  r)  =  {ze  C  |  • 

This  is  precisely  the  condition  required  for  stability  of  the  iteration, 

since  no  eigenvectors  exist  which  are  amplified  by  the  iteration.  However, 

it  is  clear  that  if  dt  is  too  large  A(L^)  will  extend  beyond  S(0,1)  and 

same  disturbances  will  be  amplified.  Thus,  explicit  time-stepping  requires 

a  time-step  control . 

The  case  of  implicit  iteration  is  more  interesting;  Equation  82  shows 

that  A(d  F)  c  LHP  =>  A(L^)  c  S(0,1)  regardless  of  the  spectral  radius  of 

d  F  or  of  dt.  This  is  because  l/(l-z),  the  spectral  function,  maps  LHP 
X 

into  S(l/2,l/2)  cS(0,l).  Hence,  the  iteration  is  stable  at  all  time- 
steps,  and  no  time-step  control  is  required  by  stability  considerations. 

That  this  argument  should  be  possible  without  assuming  d  F  to  have 

X 

constant  coefficients  is  truly  remarkable.  It  is  probably  possible  to 
extend  it  to  non-parabol'  operators  F  as  well. 


(31) 


(82) 


In  addition  to  describing  the  time-splitting,  the  remainder  of  Section 
6  will  describe  the  time  dependence  of  the  time  advance  equations  used  for 
each  physical  process  modeled  by  the  code  to  point  out  which  are  explicit 


and  which  ane  not.  Details  of  iterative  numerical  solutions  of  those 
that  are  imolicit  will  be  reserved  to  Section  7. 


For  purposes  of  reference,  we  will  think  of  the  various  stages  inter¬ 
mediate  to  the  major  time-split  processes  as  being  at  different  times,  Tne 
times  will  be  referred  to  as  t,,  t„,  t,,  t„,  t^,  and  t_.  The  time  at  the 

1  K  1  U  H  r 

initiation  of  each  cycle  is  tp  after  the  radiation  cooling,  t^^;  after 
thermal  diffusion,  t^;  after  magnetic  diffusion,  t^;  after  the  Lagranq^an 
hydrodynamics,  t,^;  and  after  the  convective  flux  which  finishes  the  cvde, 
t^,  which  is  the  same  as  t,  +  dt . 

Figure  7  shows  schematically  which  major  physical  processes  of  tne 
main  computational  cycle  advance  each  of  the  fundamental  variables  between 
those  times.  From  it,  one  can  see  that  density  is  only  changed  by  the  con¬ 
vective  transport,  while  the  internal  energy  is  changed  by  all  of  the  major 
processes  in  the  cycle. 

In  the  subsections  that  follow,  a  variable  superscripted  by  I,  R,  T, 

D,  H,  or  F  is  to  be  understood  as  having  remained  unchanged  since  t,,  t-, 
tj,  tp ,  t^,  or  tp  respectively.  Thus,  v  will  mean  the  hydrodynamic 
velocity  after  the  current  time-step's  advance  by  the  hydromagnet i c  pro¬ 
cess,  and  before  that  of  the  convective  process. 

Since  our  object  here  is  to  describe  the  time  differencing,  we  will 
use  the  continuous  vector  field  differentiation  symbols  to  write  the 
equations.  The  details  of  the  spatial  discretization  and  differencing  will 
be  discussed  in  Section  7. 

6.3  RADIATION  COOLING 

Radiation  cooling  is  the  first  physical  process  in  the  time  split  com¬ 
putational  cycle.  It  advances  the  internal  energy  from  time  t,  to  t^  by 

I  K 

P  ^  T  T  d 

e  =  e'  +  dt  fT^)^  (83) 
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Figure  7.  The  time-split  advance  of  the  fundamental  variables. 
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thus  it  is  explicitly  time  differenced.  No  time  step  control  is  used,  how¬ 
ever.  To  prevent  the  coding  from  producing  negative  energies,  the  cooling 
rate  is  limited  so  that  the  no  more  than  90  percent  of  the  internal  energy 
can  be  removed  in  one  time-step. 


6.4  THERMAL  DIFFUSION 


The  thermal  diffusion  is  implicitly  time  differenced.  The  computation 
begins  with  the  determination  of  a  local  time  centering  parameter  S.  It  is 
computed  using 


mi  n 


if  C-p  >  -K —  (3  •  '*■  7) 

T  2  Y  ^  mi  n  2  ^ 


3  =  <  2C 


T  ■  7 


(84) 


max 


^T  ^  77  ^®max  7^ 


where 


Ilk  ildt 

e 


^  3C^(dL)^ 


(85) 


is  the  thermal  Courant  number,  B  >  B  .  are  control  parameters  in  the 

max  min 


range  [0,1],  and  y  is  a  control  parameter  greater  than  1. 


The  temperature  is  then  advanced  from  tj  to  tj  by  iteratively  solving 


jJ  _  pi  dt 

'  '  '  ■  I^I 

p  c. 


3  V  •  V  +  (1  -  6)  V 


(86) 


for  T  .  Thus,  for  B  --  1.  and  B  .  =0.  the  differencing  is  fully 

max  min 


explicit  in  cells  where  Cj  <  0.75/y  and  fully  implicit  where  Cj  >  1.25/y. 
For  C-p  between  these  values,  the  differencing  is  of  mixed  type  which, 
nonetheless,  requires  iteration.  The  conductivity  tensor  is  held  fixed 


during  the  iteration,  even  though  it  may  be  a  strong  function  of  temper¬ 
ature.  The  iteration  is  carried  out  by  a  muUigrid  procedure  which  is 
described  in  detail  in  Section  10.  After  t"'”  converges  to  some  desired 
accuracy,  the  internal  energy  ismodified  explicitly  using 


6.5  magnetic  OIPFUSION 

The  third  physical  process  effected  in  the  main  cycle  is  that  of 
resistive  diffusion  of  the  magnetic  field.  This  is  the  second  of  three 
processes  in  the  code  which  are  implicitly  time  differenced.  This  compu¬ 
tation  begins  with  the  determination  of  a  local  time  centering  parameter, 
again  using  Equation  84  but  with  replaced  by  the  magnetic  diffusion 
Courant  number,  given  by 


The  magnetic  field  is  then  advanced  from  tj  to  t^  by  iteratively 
sol vi ng 


^  -  dt  [BV  X  (n^v  X  +  (1-6)  V  x  v  x  )] 


for  Thus,  for  6  =1.  and  8  •  =0.  the  differencing  is  fully 

explicit  in  cells  where  Cj,,  is  less  than  0.75/y  and  fully  implicit  where  it 
is  greater  than  1.25/y.  For  values  in  between  the  differencing  is  of  a 
mixed  type  which,  nonetheless,  requires  iteration.  Both  the  classical  and 
anomalous  resistivity  remain  fixed  during  the  iteration.  After  the  solu¬ 
tion  is  obtained  to  some  desired  accuracy,  the  joule  heating  is  applied  to 
the  fluid  using 


6.6  LAGRANGUN  HYDRODYNAMICS 


■'■'lis  is  the  other  major  physical  process  of  the  main  cycle  which  is 
implicitly  tine  differenced.  The  computation  begins  by  determining  the 
values  of  the  components  of  the  artificial  viscosity  tensor  a  from  values 
of  the  velocity  and  density  at  time  tj.  Details  of  this  are  reserved  for 
Section  7.  The  velocity,  density,  and  pressure  are  then  advanced  from  time 
tj  to  time  t^,  and  the  magnetic  field  from  time  t^  to  time  t^  by  itera¬ 
tively  solving 


-H  -I  dt  ? 
V  =  V  —  7  •  b 


S  =  -  1  )  5  +  +  a  ^ 

'  2  Y  Y  a6  a  B  aB 


a,  Be|l,2,3} 


H  I  H  -  -H 

p  =  p  -  dtp  7*v 


.  p'  .  rW 

'  3p 


-  H  I  , 

-  P  J 


-dt  -  7^) 


(91) 


Here,  3P/3pj^  is  the  square  of  the  adiabatic  sound  speed.  After  these  are 
approximately  solved  to  a  satisfactory  tolerance,  the  internal  energy  is 
advanced  to  t^  by  the  appropriate  pressure  heating  or  cooling  using 


e  -K  -j- 

p 


aB  aB 


’aS 


sH 

^aB 


3v 


H 


_ r 

3r  2 

,  3v"  a,» 


3r  • 


,  H  ^  H 

9v  3v 

t  32  ar  7 


3v 


3r 

H 

_z 

3r 


a,  Be(l,2} 


(92) 


39 


6.7  CONVECTIVE  TRANSPORT 


All  that  remains  of  the  full  model  are  the  terms  due  to  material 


derivative,  i.e., 


3(j) 

^  =  IT  "  -  'c>  •  ’‘a  ■  0 


for  <j>  =  p,  e,  V,  B.  These  are  the  same  as  Equations  71  in  Section  5, 
except  for  the  substitution  of 


+  V  •  ( pv)  =  0 


|£  +  V  .  7p  =  0 


Since  Equations  71  are  shown  in  Section  5  to  be  equivalent  to  Equations  72 
they  are  time  differenced  instead.  The  additional  p  V  •  v  of  Equation  94 
is  accounted  for  by  ignoring  p  and  using  p  in  a  single  place,  since 
Equation  91  shows  that  to  be  precisely  the  difference  between  them.  The 
time  differenced  forms  are 


/  p  p^dV  =  J  .  p^dV  +  dt  /  p  p^  (v  -  •  n  dA 

R^  R^  3R 


/  p  p^e^dV  =  /  ,  p^e^dV  +  dt  /  p  p^e^  (v  -  •  n  dA 

R^  R^  3R 


I  p  p^'v^ dV  =  /  .  p^T^ dV  +  dt  /  p  p^V^  (v  -  •  n  dA 

R^  3R  ^ 


/  p  pVdV  =  /  .  p^Tf^dV  +  dt  /  p  pV  (v  -  •  n  dA 

R^  R^  3R 


.VsT.V'j'.T'v'Cv'y.'.N’C 


where  v  is  tne  velocity  of  the  coordinate  system,  and  R  and  R  are  the 
region  at  time  to  tp  and  t^,  respectively. 
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7.  SPATIAL  DISCRETIZATION,  CENTERING,  AND  DIFFERENCING 


There  is  a  not  a  lot  of  a  detailed  numerical  analytical  nature,  such 
as  estimates  of  accuracy  or  stability,  that  requires  discussion  or  that 
will  be  illuminating.  However,  the  basis  for  the  numerical  analysis  must 
be  detailed.  This  includes  the  spatial  discretization,  the  centering  of 
the  discretized  physical  variables,  and  the  spatial  differencing  of  all 
operators  appearing  in  the  equations  which  describe  the  physical  processes 
given  in  Section  6.  The  boundary  condition  implementation  technique  and 
the  iterative  solution  procedures  for  the  implicitly  time-stepped  processes 
will  also  be  described  here. 

7.1  SPATIAL  DISCRETIZATION 

The  description  of  the  class  of  admissible  geometries  is  the  starting 
point  for  the  discretization  of  the  domain.  The  problem  domain,  a  block¬ 
like  complex,  is  the  image  of  a  block  complex  under  the  action  of  a  dif- 
feomorphism.  If  the  blocks  of  that  complex  are  discretized  in  the  most 
obvious  way,  a  discretization  of  the  block-like  complex  results  from  the 

di ffeomorphism.  So  B  ,.  is  replaced  with  the  collection  of  integer  grid- 

^  D 

points  it  contains 

8l  =  {(i,j)  e  I  1  <  i  5  It  1  5  J  5  I  =  (97) 

Since  the  location  of  the  domain  of  the  di ffeomorphi sm  is  immaterial,  each 
block  is  taken  to  have  lower-left  corner  at  (1,1).  As  the  continuous 
coordinate  system  is  described  by  the  di ffeomorphi sms 

(X^(n,  C,  t)  ,  Y^(n,  c.  t)),  1  =  1 . L  (98) 

as  in  Section  4,  then  the  discrete  coordinate  system  is  described  by 


(99) 


'  .^j  *  I  ^  f  M  h  M  J  ^  ^  1 . LI 

a  finite  collection,  indexed  by  1,  of  two  dimensional  arrays  of  points  in 
the  computational  plane.  These  arrays  of  points  are  referred  to  as  the 
grid,  or  the  mesh,  in  which  context  they  are  thought  of  as  being  connected 
to  their  nearest  neighbors  which  have  the  same  value  of  i,  or  of  j,  as 
shown  in  Figure  3.  The  four  points 


(''Ij-’lji  (*l  *  1.  j- ’1  *  1,  j) 

j  +  l’  ''i.  j  ♦  l'  ("I  ♦  1,  j  ♦  ’i  *  1.  j  *  l' 


(100) 


which  lie  in  general  position,  and  the  four  lines  of  constant  i  or  constant 
j  which  connect  them  form  the  cell  indexed  by  The  location 


’’-j)  *  (*!  ♦  i.j  t'  *  i.j) 

*’  f  *1  ,j  ♦  r  ''i  ,j  *  1^ 


^^i  +  1,  j  +  1‘  ^  +  1,  j  + 


(101) 


is  called  the  cell  center,  and  points  of  the  arrays  are  known  as  vertices. 
Such  a  mesh  is  often  referred  to  as  an  arbitrary  quadrilateral  mesh,  for 
obvious  reasons. 


All  of  the  variables  in  the  code  which  have  spatial  dependence  are 
stored  as  functions  of  1  and  (i.j),  so  that  the  density,  for  example,  is 
known  as  a  function  of  (x,y)  only  through  the  parametric  description  given 
by  Equation  99  and  the  collection  of  arrays 

loljl  (102) 


44 


7.2  SPATIAL  CENTERING 


For  the  purpose  of  deriving  difference  formulas,  most  spatially  vary¬ 
ing  quantities  are  thought  of  as  having  values  located  at  either  cell 
centers  or  vertices.  The  fundamental  variables  of  the  code  are  thus  split 
into  two  classes; 

Vertex  Centered  --  position,  velocity 

Cell  Centered  --  density,  internal  energy,  magnetic  field 

There  are  many  other  space  dependent  ouantities  derived  from  these,  and 
they  are  usually  quite  easily  identit  ;  as  belonging  to  one  or  the  other 
class  by  the  following  simple  rule:  difference  operators  map  each  class  to 
the  other.  Thus,  the  current  density  is  a  vertex  centered  quantity,  since 
0=7x8".  When  it  is  necessary  to  call  attention  to  their  centering,  cell 
quantities  will  be  indicated  by  a  superscript  "c",  and  vertex  quantities  by 
a  superscript  "v". 

Unfortunately,  a  third  class  of  variable  exists:  the  edge  centered 
quantity.  The  convective  fluxes  of  mass,  momentum,  internal  energy,  and 
magnetic  field  fall  into  this  class;  they  are  all  derived  quantities.  When 
it  necessary  to  call  attention  to  such  a  quantity's  centering  it  will  be 
superscripted  with  an  "e". 

7.3  FINITE  VOLUME  DIFFERENCING 

In  general,  the  finite  volume  approach  is  used  to  derive  the  differ¬ 
encing  in  the  code.  This  approach  can  best  be  understood  by  considering 
the  example  of  7  x  (  ).  One  form  of  Stokes  Theorem  that  involves  this 
operator  is 


fr,  7  X  F  dV  =  n  X  F  dA 


fl03) 


Here  R  is  region  of  3-space  and  3R  is  its  bounding  surface.  If  R  is  taken 
to  be  the  small,  discrete  volume  in  3-space  described  by  a  single  cell  in 
the  computational  grid,  then  Equation  103  becomes 

[7  X  f]\  dvj  .  =  y  (n  X  F)^  dAj  (104) 

^  ^  faces  ij  ^ 

where  it  has  been  assumed  that  v  x  F  and  n  x  F  are  uniform  over  the  volume 

and  its  faces,  respectively.  The  quantity  dv|.  here  is  the  cell  volume 

and  dA  ,  the  face  area.  If  values  for  n  x  F  are  specified  for  each  of 
i  j  _  _ 

the  faces  of  the  cell,  then  Equation  104  can  be  used  to  define  v  x  F.  The 

differences  are  hidden  in  the  direction  of  the  normals  in  the  sum  on  the 

right  hand  side;  the  normal  on  one  face  is  opposed  in  sense  by  the  normal 

on  the  opposite  face.  In  what  follows,  the  indices  [  will  be  omitted. 

i  j 


The  principal  advantage  of  this  scheme  of  differencing  is  that  conser¬ 
vation  laws  written  in  terms  of  the  quantities  in  the  vector  integral  theo¬ 
rems  are  usually  well  respected  by  the  differencing.  In  fact,  for  the 
example  above,  all  that  is  required  to  insure  that  Equation  104  is  satis¬ 
fied  exactly  for  R  a  region  composed  entirely  of  whole  cells  is  that  the 
quantity  rT  x  F  be  computed  the  same  way,  except  for  its  sign,  for  a  par¬ 
ticular  face  regardless  of  the  cell  with  which  it  is  associated.  This  is 
worthy  of  more  explication.  First  note  that 


I  (7  X  F)  dV  =  I 
cells  cells 


I  (n  X  F)  dA 
faces 


;io5) 


Since,  by  the  above  assumption,  two  terms  (n  x  F)  dA  appear  within  the 
double  sum  for  each  cell  interface,  each  with  opposing  sign,  all  but  the 
contributions  from  the  exterior  faces  must  cancel.  Hence  it  is  true, 
exactly,  without  approximation,  that 


I  (V  X  F)  dV  =  I  (n  X  F)  dA 
cells  exterior 

faces 
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This  is  the  best  possible  statement  of  Equation  103  for  the  discretized 
data  available. 


There  are  two  subtleties  in  the  application  of  the  finite  volume  dif¬ 
ferencing  to  any  particular  case,  however.  The  first  is  due  to  the  pos¬ 
sible  cylindrical  symmetry.  The  second  is  due  to  the  different  centering 
types  of  the  data.  There  is  also  a  bit  of  unpleasantry:  the  requirement 
to  specify  (n  x  F)  dA  precisely  in  each  case.  Each  of  these  shall  be  dealt 
with  in  turn. 

In  cylindrical  symmetry,  though  all  vector  components  satisfy 

^  F,  =  0  (107) 

the  standard  basis  (r,  e,  z)  satisfies 


ii  -  -  - 
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If  the  three-dimensional  discrete  volume  associated  with  a  two-dimensional 
cell  in  (r,z)  space  is  assumed  to  extend  a  small  increment  de  in  the  9 
direction,  as  shown  in  Figure  9,  then  Equation  104  becomes 
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Figure  9.  Finite  volume  for  cell  centered  differences  of  vertex  quantities 
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since  F  is  independent  of  9.  Here  da  is  the  area  of  the  cell  in  the  com 
putational  plane,  and  the  sum  extends  over  the  four  faces  obtained  by 
extending  the  cell  edges  normal  to  the  computational  plane  by  the  incremen 
rde  .  Taking  the  limit  as  de  goes  to  zero  and  using  Equation  108  we  obtain 
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Thus  the  difference  formula  for  7  x  F  contains  terms  in  cylindrical  sym¬ 
metry  that  are  absent  in  planar  symmetry.  In  the  code  planar  symmetry  is 
easily  implemented  by  setting  the  quantity  da  to  zero  and  the  radius  r  in 
the  summation  in  Equation  110  to  one. 

It  might  seem  that  all  occurrences  of  the  curl  operator  could  be  dif¬ 
ferenced  alike;  sadly,  this  is  not  the  case.  There  are  two  distinct  dif¬ 
ference  forms  for  it  in  use  in  the  code: 

cell  centered  curl  of  vertex  quantity 
vertex  curl  of  a  cell  centered  quantity 


and  this  problem  extends  to  other  operators  as  well.  In  large  measure  the 
cases  differ  only  in  the  choice  of  the  volume  over  which  the  technique  is 
applied.  The  finite  volume  in  Figure  9  is  the  one  appropriate  for  comput¬ 
ing  cell  centered  differences  of  vertex  centered  quantities.  Figure  10 
shows  a  region  in  the  computational  plane  which,  when  extended  by  rde  in 
the  0  direction,  becomes  the  finite  volume  that  is  used  to  compute  vertex 
centered  differences  of  cell  centered  quantities.  The  region  consists  of 
four  parallelograms  formed  by  the  vertex  (i,j)  and  the  midpoints  of 
adjacent  pairs  of  edges  leaving  it.  Since  this  is  the  most  common  form  of 
differencing  in  the  code,  some  geometric  quantities  associated  with  this 
finite  volume  are  computed  once  per  cycle  and  saved  for  reuse  to  avoid 
recomputation. 


7.4  STORED  GEOMETRIC  COEFFICIENTS 

Those  geometric  quantities  are  most  easily  understood  in  the  context 
of  the  derivation  of  the  vertex  centered  gradient  of  some  cell  centered 
scalar  quantity  > .  Applying 
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where  n  is  an  outward  pointing  vector  field,  to  the  finite  volume  described 
above  yields 
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Here  the  first  sum  is  over  the  four  shaded  paral 1  el ograms  and  the  second  is 
over  the  eight  edges  of  the  shaded  figure  of  Figure  10. 

The  second  sum  is  collapsed  to  four  terms  and  converted  to  a  sum  over 
adjacent  cells  by  assuming  that  within  each  cell  4.  is  constant.  The  two 
terms 


f  dA>  f  dA> 


(113) 


within  the  sum  which  derive  from  a  single  cell,  shown  in  Figure  11,  are 
replaced  with  the  single  term 


^n  4^1 


(llA) 


where  the  subscripts  t,  s,  and  d  stand  for  top,  side,  and  diagonal  as 
illustrated  in  the  figure.  This  is  equivalent  to  assuming 


which  is  not  precisely  true  in  cylindrical  symmetry,  since  it  ignores  the 
variation  in  dA  with  r  within  the  cell.  This  does  not  lead  to  disastrous 
results  since  we  are  using  the  finite  volume  approach.  The  boundary 
quantities  need  only  be  estimated,  not  computed  exactly,  since,  so  long  as 
they  are  all  computed  consistently,  the  integral  conservation  laws  remain 
valid.  The  same  radius  is  used  for  each  of  the  four  directed  vertex  face 
areas  in  one  cell;  that  radius  is  the  centroid  of  the  cell  area.  These 
areas  are  thus  given  by 
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where  r^  represents  the  cell  centroid,  R  is  a  rotation  by  v/2  radians,  dl 
is  the  diagonal  of  the  corner  parallelogram  as  indicated  in  Figure  11  and 
the  corner  index  k  is  as  shown  in  Figures  8  and  11. 

The  four  vector  quantities  in  Equation  116,  and  the  quantity 
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referred  to  as  the  vertex  volume,  and  the  areas  da^  for  k  =  1 . 4  of 

parallelograms  such  as  the  one  in  Figure  11,  completely  specify  the  vertex 
differences  of  cell  centered  quantities.  They  are  computed  once  per  cycle 


and  stored  for  reuse.  The  corner  areas  da.  and  the  directed  vertex  areas 


dA  c 


(n  are  cell  quantities,  as  indicated  by  the  superscript  c,  but  are 
each  associated  with  a  particular  vertex  of  the  cell,  and  are  numbered  as 
the  vertices  are  numbered  in  Figure  8.  In  addition  to  these,  the  four 
volumes 
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referred  to  as  the  corner  volumes,  and  the  reciprocal  cell  volume 
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are  stored  for  use  in  other  differencing  and  averaging.  Note  that  the 
vertex  volume  is  defined  as 


(f)'-  J 

adjacent 

cells 

where  the  notation  k(c)  indicates  that  the  appropriate  corner  number  is  to 
be  chosen  within  each  cell  so  that  the  corner  referred  to  is  adjacent  to 
the  vertex  at  which  the  volume  is  being  computed.  One  advantage  of  comput¬ 
ing  these  geometric  quantities  in  only  one  place  in  the  code  is  that  they 
can  all  be  easily  changed  in  a  consistent  way  if  that  should  be  necessary. 
That  has  been  done  on  at  least  one  occasion. 


With  the  above  definitions,  the  difference  form  of  the  vertex  centered 
gradient  of  a  cell  centered  scalar  quantity  ()>  may  be  written  as 


[- 


.  I 

I  adjacent 
(  cells 


4)  [-r  da 


k(c) 


(121) 


55/56 


%  s 


8.  DIFFERENCE  OPERATORS  BY  PHYSICAL  PROCESS 


The  remainder  of  this  section  will  follow  the  organization  of  the  main 
computational  cycle,  dealing  with  each  type  of  differencing  for  each  oper¬ 
ator  where  it  first  makes  an  appearance.  The  precise  specification  of  the 
boundary  quantities  for  non-vertex-centered  difference  operators  will  also 
be  covered  along  the  way,  since  they  depend  on  the  centering  type  of  the 
data  to  which  the  operator  applies. 

8.1  RADIATION  COOLING 

There  is  no  differential  operator  in  Equation  83  which  governs  the 
radiation  cooling  model.  That  is,  spatial  derivatives  play  no  role  in 
energy  loss  by  free-free  emission. 

8.2  THERMAL  DIFFUSION 

There  are  two  differential  operators,  the  divergence  and  the  gradient, 
in  the  thermal  diffusion  model  of  Equation  86.  The  divergence  must  produce 
a  cell  centered  quantity,  the  thermal  heating  rate,  and  the  gradient  is  of 
the  temperature,  a  cell  centered  quantity.  The  centering  of  the  gradient 
is  arbitrary,  so  it  may  be  taken  to  be  a  vertex  quantity. 

Since  the  gradient  of  the  temperature  is  a  vertex  difference  of  a  cell 
centered  quantity,  it  is  computed  using  Equation  121  with  ij>  =  T  and  the 
geometric  quantities  defined  above.  The  thermal  flux,  also  a  vertex 
quantity,  is  computed  by  multiplying  by  the  effective  thermal  conductivity, 
so  that 


r  =  k  7  T 


(122) 
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The  divergence  of  the  flux  is  then  computed  to  determine  the  cell 
centered  thermal  heating  rate.  Because  the  thermal  flux  is  a  vertex 
centered  quantity  and  the  thermal  heating  rate  is  cell  centered,  the  geo¬ 
metric  quantities  do  not  apply  exactly  as  derived.  However,  the  boundary 
integral  areas  for  vertex  centered  differences  are  related  to  those  for 
cell  centered  differences.  This  relationship  is  made  clear  in  Figure  12. 
The  assumption  there  is  that  the  areas  of  the  surfaces  associated  with  the 
cell  edges  do  not  vary  with  radius.  To  correct  for  this  error,  additional 
terms  proportional  to  the  cell  corner  areas  are  included  in  the  dif¬ 

ference  formula.  Thus,  the  divergence  of  thermal  flux  is  computed  using 
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8.3  MAGNETIC  DIFFUSION 

The  only  differential  operator  in  Equation  89,  which  governs  diffusion 
of  magnetic  field,  is  the  curl.  However,  two  centering  types  are  required. 
As  mentioned  above,  since  B"  is  cell  centered,  J  is  taken  to  be  vertex 
centered.  Hence,  the  difference  operator  in 


J  =  V  X  B 


(124) 


is  a  vertex  centered  difference  of  a  cell  centered  quantity,  while  that  in 
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is  a  cell  centered  difference  of  vertex  data. 

The  derivation  of  the  difference  formula  for  the  cell  centered  curl  of 
a  vertex  centered  quantity  was  begun  above,  in  the  introductory  discussion 
on  finite  volume  differencing.  All  that  remains  to  complete  the  derivation 
is  to  specify  the  quantities  on  the  right  hand  side  of  Equation  110, 
repeated  here  with  F  =  n  •  J  substitued  for  F 
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for  easy  reference.  The  volume  on  the  right  hand  side  needs  to  be  the 
volume  of  the  cell  per  radian,  and  hence  the  reciprocal  cell  volume,  one  of 
the  saved  geometric  quantities,  is  used  there.  Since  E  is  vertex  centered, 
the  first  term  in  the  braces  must  be  expanded  into  an  average  over  the  four 
vertices  of  the  cel  1 
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The  da^  here  are  the  cell  corner  areas,  which  are  also  among  the  previously 

defined  geometric  quantities,  and  the  notation  v(k)  indicates  that  the 

^  ,•  indices  of  E  are  chosen  to  match  the  corner  index  k  as  in  Figure  8. 
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The  terms  of  the  sum  which  remains  are  first  rearranged  as 
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The  first  factor  is  defined  by 


(128) 
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where,  as  before,  the  primed  and  unprimed  quantities  refer  to  the  two 
vertices  which  are  the  ends  of  the  appropriate  edge.  The  second  factor  is 


taken  to  be 

^"ndl  =  RdT  ( 13C 

in  which  R  is  the  r/ 2  rotation  operator,  and  dl  is  the  edge,  oriented  so 
that  the  cell  boundary  is  traversed  counterclockwise. 
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The  vertex  centered  curl  of  cell  centered  data  is  derived  using  the 
finite  volume  described  by  Figure  10.  It  is  given  by 
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where  all  the  geometric  quantities  have  already  been  defined. 

8.4  LAGRANGIAN  HYDRODYNAMICS 

Here  the  differential  operators  in  the  equations  of  advance  are  the 
divergence  and  tne  covariant  derivative,  that  is 


V  .  (  )  and  (  )  •  V  (  )  (132 

The  divergence  is  required  in  two  centering  types,  while  the  covariant 
derivative  must  produce  a  cell  centered  result  from  vertex  centered  data. 

First,  the  momentum  equation  requires  a  vertex  centered  divergence  of 
a  cell  centered  tensor  of  rank  two.  This  difference  formula  can  thus  be 
written  in  terms  of  the  saved  geometric  quantities.  The  form  is 
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which  differs  from  Equation  131  by  the  substitution  of  the  cross  product 
for  the  dot  product,  and  since  F  is  a  tensor  instead  of  a  vector. 

The  second  occurrence  of  the  divergence  here  is  in  the  terms  repre¬ 
senting  the  effects  of  the  velocity  field  on  the  magnetic  field  and  mass 
density.  Both  of  these  require  a  cell  centered  divergence  of  a  vertex 
centered  quantity,  the  velocity.  The  divergence  of  the  velocity  is  done 
precisely  the  same  as  the  divergence  in  the  thermal  transport,  that  is. 
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The  covariant  derivative  required  is  S  ♦  Tv,  and  by  a  similar  argument 
to  that  given  above,  is  given  by 
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R.5  CONVECTIVE  TRANSPORT 


The  integral  statements  of  the  convective  process.  Equations  96  of 
Section  6,  can  be  differenced  directly.  Application  of  the  first  of  those 
integrals,  which  governs  the  flux  of  mass,  to  a  single  cell  produces 
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The  fractional  timesteps  of  the  cell  volumes  have  to  be  displayed  since  the 
grid  movement  step  is  combined  with  the  convective  step.  No  terms  due  to  w 
appear  here  since  the  density  is  independent  of  9.  When  the  convection  of 
the  vector  quantities  8"  and  v  is  considered  below,  the  turning  of  the  basis 
vectors  will  require  treatment  of  such  terms.  The  edge  quantities 
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which  are  the  rate  of  flux  of  fluid  volume  across  the  edges,  are  taken  to 
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and 
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The  speci f ication  of  the  edge  value  of  the  density  determines  the 
order  of  accuracy  and  the  stability  of  the  difference  scheme  for  its  trans¬ 
port.  In  general,  techniques  for  selecting  the  edge  value  are  based  on 
donor  cell  averaging,  which  is  of  first  order  accuracy,  and  linear  interpo¬ 
lation,  which  is  of  second  order.  The  order  of  accuracy  referred  to  here 
is  that  relevant  to  smooth  functions,  i.e.,  those  for  which  the  relative 
change  from  cell  to  cell  is  small.  In  such  a  circumstance  second  order 
differencing  is  clearly  superior.  However,  if  the  relative  cell  to  cell 
change  is  large,  both  schemes  produce  errors.  Second  order  differencing 
tends  to  introduce  unphysical  behavior  in  the  solution  while  first  order 
differencing  tends  to  spread  short  scale  length  phenomena.  Since  neither 
scheme  can  support  scale  lengths  shorter  than  the  mesh  size,  it  is  usually 
felt  that  first  order  differencing  is  superior  in  these  regions.  Thus 
these  schemes  are  often  combined  in  some  way,  in  an  effort  to  increase 
accuracy  where  the  solution  is  smooth  and  decrease  unphysical  behavior 
where  it  is  not. 


Flux  corrected  transport,  described  by  Zalesak  (Ref.  9),  is  one 
approach  to  this  averaging.  Another  approach  is  to  form  a  combination  of 
the  two  schemes,  in  effect,  an  average  of  the  averages.  If  the  weights  are 
chosen  adaptively,  as  some  function  of  the  computed  solution,  the  combined 
difference  scheme  can  be  made  second  order  accurate  where  the  solution  is 
smooth  and  first  order  where  it  is  not.  This  is  the  approach  followed  by 
Brackbill  in  MOQUI,  and  adopted  here.  Letting  the  superscripts  i  and  o 
stand  for  inside  and  outside,  the  edge  value  of  the  density  is  given  by 
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where  the  function  a  that  determines  the  order  of  the  differencing  is 
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Here,  V  is  the  edge  volume,  dV  is  the  edge  volume  flux  given  in  Equation 
137,  and  s  =  sign  (dV  ).  Where  3=1,  the  edge  density  p  is  determined  by 
linear  interpolation  between  p’  and  p°.  Where  3=0,  it  is  either  p^  or 
depending  only  on  whether  the  volume  flux  is  into  or  out  of  the  cell  in 
question,  and  hence  the  differencing  will  be  donor  cell.  For  values  of  3 
in  (0,1)  the  difference  scheme  is  a  linear  combination  of  the  two.  The 
control  parameter  s  is  determined  by 
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The  function 


3^  =  1  -  ^  (h(g^  )  +  h(g°))  (143) 

is  a  measure  of  the  local  gradient  of  the  relative  velocity  of  the  grid  and 
the  fluid,  with  h  as  shown  in  Figure  12,  and  g  computed  as 
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where  d  is  measure  of  the  local  grid  spacing.  The  choice 
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for  the  parameters  of  Figure  13  would  cause  donor  cell  differencing  to  be 
used  wherever  the  vertex  to  vertex  variation  of  the  relative  velocity  of 
the  grid  and  fluid  exceeds  10  percent.  Of  course,  the  choice 
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Figure  13.  Convective  flux  differencing  control  function. 
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forces  donor  cell  differencing  to  be  used  everywhere.  This  is  the  default 
choice,  and  the  only  one  which  has  ever  been  successfully  used  for  real 
s i nu 1  at i ons  . 


■^he  other  control  function  in  Equation  142  is  a  measure  of  the  density 
gradient.  It  is  computed  as 
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where  s^  =  sign(p°/p^  -  1  ).  The  reader  is  invited  to  show  that  this  pro¬ 
duces  1  if  the  velocity  of  the  grid  relative  to  the  fluid  is  toward  larger 
density,  and  p ^iha  1 1 large  -jf  velocity  is  toward  smaller  density.  The 
latter  value  has  the  effect  of  limiting  the  edge  density  to  (3/2) 
when  the  ratio  of  the  densities  is  very  small  and  the  grid  is  moving  into 
the  lower  density  fluid,  rather  than  the  allowing  the  large  value  that 
linear  interpolation  between  and  would  produce. 


8.5.1  Homogeneity  of  the  Transport  Scheme — It  is  desirable  for  the 
mass  transport  process  to  have  the  following  property; 

Homogeneity;  For  0  fluid  velocity  and  uniform  density,  and  for  any 
grid  velocity,  the  density  after  the  transport  step  should  be  the  same 
as  the  density  before. 


Letting  p  be  the  value  of  the  uniform  density,  it  is  clear  from  Equation 
^He 

140  that  p  =  0^.  Equations  136  and  137  then  imply  that 
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Hence  o 


if  and  only  if 
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that  is,  if  only  if  the  estimates  of  volume  flux  due  to  coordinate  system 
notion  exactly  account  for  the  changes  in  cell  volume  due  to  grid  motion. 
Unfortunately,  the  volume  flux  estimates  given  by  the  expressions  above  do 
not  satisfy  this.  The  errors  are  of  two  kinds.  First,  the  fluxes  of 
volume  to  the  diagonal  neighbor  cells,  are  not  included  in  the  estimates. 
Second,  even  when  the  grid  moves  so  that  the  diagonal  fluxes  are  zero,  the 
edge  volume  flux  estimates  are  in  error  due  to  incorrect  radius  weighting 
and  improper  time  centering.  The  complexity  of  making  the  volume  fluxes 
exact  is  deemed  to  be  a  greater  penalty  than  the  compensation  that  would 
result . 

a. 5. 2  Internal  Energy  Transport--The  other  transport  equations  of 
Equation  96  can  be  made  homogeneous  by  the  use  of  mass  flux  based  trans¬ 
port.  For  example,  consider  the  spatially  discretized  version  of  the 
second  of  Equations  96,  which  governs  energy  transport, 
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If  the  cell  mass  m  is  defined  as 


c  rdV^,' 


and  the  edge  mass  flux  dm  as 


.  G  Hg  ( —  >  F  dA  ..  Hg 

dm  =  p  [v^  -  V  J  •  n  dt  =  p  dv 


then  Equation  150  can  be  written  as 


F  Fc  H  Ic  ^ 
e  m  =  e  m  + 


r  e  .  e 
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u 

edges 


These  edge  mass  fluxes  are  precisely  the  same  as  those  computed  during  the 
mass  transport  phase.  To  emphasize  this  fact  note  that  Equation  136  can  be 
wri tten 


Fc  Ic 
m  =  m 


I 

edges 


(154) 


Thus  the  mass  fluxes  computed  there  are  saved  and  reused  during  the  trans¬ 
port  of  all  other  quantities.  The  edge  value  of  the  internal  energy  is 
computed  by  the  same  scheme  as  that  for  the  density.  In  fact,  the  calcu¬ 
lation  in  Equations  140  through  147  is  repeated  exactly  with  the  exception 
of  the  substitution  of  e  for  p  throughout.  Now,  if  e^  =  e^,  then,  just  as 
for  p  above,  e®=  e°.  In  this  case  Equation  153  becomes 


e=e„m  +  Y  dm  1/m 

0  1  4-  i 

edges 


which,  by  Equation  154,  implies  e  =  e^. 


8.5.3  Monotonicity  of  the  Transport  Scheme--Another  significant 
property  for  a  transport  scheme  is  monotonicity.  A  scheme  is  monotonic  if 
the  transport  of  a  quantity  from  one  cell  to  a  neighbor,  in  the  absence  of 
any  other  transport,  results  in  the  final  densities  of  that  quantity  in 
those  cells  lying  between  the  initial  densities,  and  in  the  original  order. 
The  mass  flux  based  transport  scheme  described  above  is  monotonic  if  donor 
cell  differencing  is  used,  though  distinctly  non-monotonic  if  any  linear 
interpolation  is  used.  This  will  demonstrated  for  the  transport  of 
internal  energy.  The  superscripts  +  and  -  be  used  to  designate  quantities 
belonging  to  the  cells  to  which  mass  is  added  and  from  which  mass  is 


removed,  respectively.  With  this  notation,  donor  cell  differencing  corre- 

e  H- 

sponds  to  the  choice  e  =  e  .  For  the  case  of  mass  transport  toward 
higher  energy,  the  monotonicity  property  is  easily  stated  as 


H-  H  + 

e  <  e 


<  e^-  <  e*"^  <  e^^ 


Since  there  is  mass  flux  across  only  one  edge  Equation  155  becomes 
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Because 


H-  ,  e  ,  H+ 
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(158) 


it  is  clear  that 


e«-  <  d"' 

m^ 


(159) 


p  H- 

with  equality  on  the  right  only  if  e  -  e  .  The  remainder  of  the  ine¬ 
qualities  required  to  prove  monotonicity  for  donor  cell  differencing  follow 
in  an  analogous  fashion.  That  strict  inequality  on  the  left  in  Equation 
158  implies  strict  inequality  on  the  right  of  Equation  159  shows  that  non- 
donor  cell  differencing  produces  non-monotonicity.  When  a  transport  scheme 
is  non-monotonic,  a  propagating  front  separating  two  ostensibly  uniform 
states  will  have  excursions  above  and  or  below  the  limit  values.  Such  an 
excursion  is  particularly  troublesome  when  it  causes  a  sign  change  in  the 
energy!  It  is  for  this  reason  that  donor  cell  differencing  is  used  almost 
exclusively  on  real  problems. 

8.5.4  Cylindrical  Effects  from  Convective  Derivati ves--The  last  two 

integrals  of  Equation  96,  having  vector  integrands,  generate  terms  pro- 

portional  to  w  -  w  .  Of  course,  the  coordinate  system  does  not  move  in  the 
*  c  l_l 

6  direction,  i.e.,  w^  =  0,  and  thus  these  terms  are  proportional  to  w  . 

Consider  v  •  7F  in  cylindrical  symmetry 


w 

'w 
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Since  the  coordinate  functions  F^,  F^ ,  and  F^  are  independent  of  9  by  the 
symmetry  assumption,  only  the  basis  vectors  contribute  to  the  9  derivative 
term  and  thus 


■■  "F  =  (u  ^  +  v  F  +  -  (F  9  -  F  r) 
3r  az-*  r  r  9 


For  F  =  V,  these  terms  are  responsible  for  conservation  of  angular  momentum 
and  centripetal  acceleration,  while  for  T  =  B"  ,  they  account  for  the  turn¬ 
ing  of  trapped  field  lines  due  to  fluid  rotation.  They  seem  more  appropri¬ 
ately  included  with  the  Lagrangian  hydrodynamics,  and  appear  there  in 
explicit  form.  These  terms  are  omitted  from  the  difference  formulas  in  the 
convective  processes. 

8.5.5  Magnetic  Flux  Transport--The  difference  form  for  the  fourth  of 
Equations  96,  which  governs  the  transport  of  magnetic  field,  is  thus 
written 


=  fm^''  B^  +  Y  B^  dm®}  /  m^"’' 
edges 


(162) 


Remember  that  the  cell  mass  appears  here  to  remove  the  effect  of  the 
divergence  of  the  fluid  velocity  which  is  properly  included  in  the  differ¬ 
ence  equations  of  the  Lagrangian  hydrodynamics.  The  edge  value  ^  is  given 
by  Equations  140  and  141  with  p  replaced  by  B”  and  b  =  6^  of  Equation  143. 

8.5.6  Momentum  Transport — The  transport  of  momentum,  described  by  the 
third  of  Equations  96,  differs  from  the  other  three  transport  processes 
described  above  because  the  velocities,  which  play  the  role  of  densities, 
are  vertex  rather  than  cell  quantities.  One  possibility  is  to  form  ind 


transport  the  cell  centered  momentum  using  the  ceil  mass  and  the  velocity 
averaged  from  the  four  vertices,  and  then  reform  the  vertex  velocity  after 
the  transport  step.  That  is  very  diffusive,  since  it  transfers  momentum 
across  a  cell  even  in  the  absence  of  mass  flux.  Brackbill's  approach  in 
MOQUI  was  to  use  a  separate  vertex  control  volume  similar  to  that  of  Figure 
8,  and  to  transfer  momentum  on  the  basis  of  the  mass  fluxes  between  those 
volumes.  Those  difference  equations  suffered  from  lack  of  the  homogeneity, 
and  hence  a  grid  moving  through  a  uniform  velocity  field  introduced  non¬ 
uniformity.  The  problem  was  caused  by  the  different  mass  fluxes  used  for 
mass  and  momentum  transport. 


It  is  possible  to  avoid  both  of  these  problems  at  the  expense  of  some 
increased  computational  effort  and  additional  complexity.  The  scheme 
req'uired  is  baseo  on  the  transport  of  four  cell /vertex  momenta  for  each 
component  of  the  velocity.  For  each  vertex  indexed  by  k  =  1,...,4  those 
momenta  are  given  by 
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(163) 


The  flux  step 


edges 


(164) 


uses  exactly  the  same  mass  fluxes  as  for  the  transport  of  mass, 
flux  step  is  performed  the  new  vertex  mass  is  computed  from 


After  the 
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adjacent 
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Then  the  new  vertex  momentum 


^  1 


J  ^  ^(c) 

adjacent 


(166) 


is  obtained,  and  from  it  and  the  vertex  mass,  the  new  vertex  velocity  is 
obtained  by 


— Fv  ,  vF  / ,  r-7  \ 

V  =  P  /m  (167) 

In  the  cylindrical  case,  the  9  velocity  is  first  replaced  by 

-  u^dtjw^*^  (168) 


and  divided  by  r^  after  Equation  167  is  performed. 

This  scheme  can  easily  be  shown  to  conserve  total  momentum  exactly. 
The  precise  statement  which  can  be  proven  is 


I  =  I  ,  I  ^  1  Tj)®  dm^  (169) 

all  all  exterior  (c=l 

vertices  vertices  edges 

That  is,  the  change  in  total  momentum  is  caused  solely  by  fluxes  through 
the  problem  boundaries.  Furthermore,  when  there  is  no  mass  flux,  there  is 
no  change  in  the  velocity. 

8.5.7  Energy  Conservation  of  the  Transport  Scheme — The  total  energy 
is  not  conserved  by  the  transport  scheme  described  above.  The  error  occurs 
because  neither  magnetic  energy  nor  kinetic  energy  is  conserved  exactly. 

The  change  in  kinetic  energy  due  to  the  transport  scheme  can  be  estimated 
by  considering  the  transport  of  mass  and  momentum  due  to  the  relative 
motion  of  the  grid  and  fluid  at  a  single  vertex.  Figure  14  shows  a  uniform 
rectangular  grid  near  a  central  vertex,  and  a  new  location  for  that  vertex 
and  its  associated  edges.  The  mass  fluxes,  dm,  which  this  grid  motion 
causes  are  also  shown;  the  cell  mass  m  is  the  same  in  all  cells.  The 
velocity  v  is  assumed  directed  vertically  downward  at  all  vertices.  This 


velocity  is  taken  to  be  constant  along  the  horizontal  grid  lines,  but  a 
constant  ratio  y  is  assuned  between  velocities  at  vertices  connected  by  a 
vertical  edge.  The  ratio  is  taken  in  such  a  way  that  if  the  velocity 
decreases  in  the  direction  of  the  mass  flux  shown  in  the  figure,  then 
Y  >  1.  The  motion  of  this  single  vertex  changes  the  momentum  of  all 
9  vertices  shown  in  the  figure.  Assuming  donor  cell  averaging,  the  rela¬ 
tive  change  in  total  kinetic  energy,  normalized  to  9  times  the  kinetic 
energy  of  the  central  vertex,  is 

(y  - 1)^  [1^-2)  (4  "  k  ^  0  ^"4}  (170) 

Y 

where  e  =  dm/m.  Thus  kinetic  energy  is  lost  where  the  grid  moves  toward 
higher  velocity,  and  gained  where  it  moves  toward  lower  velocity. 

It  is  possible  to  force  the  conservation  of  total  energ[y  by  transport¬ 
ing  it  instead  of  internal  energy,  as  described  by  Brackbill  in  (Ref,  5). 
This  scheme  is  optional  in  the  code.  When  it  is  used,  a  loss  of  kinetic 
energy  during  the  momentum  and  mass  transport  steps  will  result  in  a 
decrease  in  the  internal  energy.  The  integrated  effect  may  be  large  enough 
to  cause  the  internal  energy  of  some  cell  to  become  negative.  In  addition 
to  being  unphysical,  this  results  in  interpolation  off  the  edge  of  the 
equation  of  state  tables,  and  sometimes  produces  negative  temperature  or 
negative  squared  sound  speed.  As  the  square  roots  of  these  quantities 
figure  in  the  computation  of  Spitzer  di f fusi vi ties  and  time-step, 
respectively,  this  will  cause  the  simulation  to  terminate  with  a  floating 
point  error. 

The  principal  cause  of  these  problems  is  that  the  mass  fluxes  and  edge 
densities  of  the  other  fluxed  quantities  are  selected  explicitly.  A  trans¬ 
port  scheme  which  conserves  total  energy,  momentum,  and  mass  and  maintains 
positivity  of  internal  energy  and  density  probably  requires  the  simultane¬ 
ous  implicit  determination  of  the  fluxes  of  these  quantities  from  the 
differenced  forms  of  the  conservation  laws.  Such  a  procedure  would  be 
expensive  to  code  and  to  run. 
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9.  DISCRETIZATION  OF  BOUNDARY  CONDITIONS 


Each  difference  equation  detailed  in  the  previous  sections  applies  at 
any  vertex  or  cell  where  all  data  its  evaluation  requires  is  defined  at  all 
surrounding  vertices  or  cells.  Clearly,  there  niust  be  edges  or  boundaries 
since  only  a  finite  amount  of  data  may  be  stored.  Thus  some  data  around 
the  edge  of  the  problem  region  must  be  generated  by  separate  equations; 
these  are  the  difference  forms  of  the  physical  boundary  conditions.  There 
are  two  distinctly  different  ways  they  could  be  applied. 

In  the  first,  the  physical  boundary  conditions  are  applied  to  deter¬ 
mine  expressions  for  the  data  outside  but  adjacent  to  the  edge  of  the 
problem  region  in  terms  of  data  inside  the  problem  region.  Then,  different 
difference  equations  for  the  edges  are  derived  by  substituting  these 
expressions  into  the  general  difference  equations. 

The  second  approach  avoids  the  derivation  of  the  boundary  difference 
equations;  space  is  provided  for  the  surrounding  data,  and  those  values  are 
computed  explicitly  from  the  expressions  derived  from  the  physical  boundary 
conditions.  Thus  the  full  difference  equations  may  be  applied  to  the 
boundary  cells  or  vertices,  just  as  to  the  interior  cells  and  vertices. 

The  extra  storage  locations  are  referred  to  as  ghost  cells  or  ghost 
vertices,  depending  on  the  type  of  the  data  under  consideration,  and  this 
approach  is  called  the  ghost  cell  technique.  The  cells  and  vertices  where 
the  full  difference  equations  are  applied  are  referred  to  as  real  cells  or 
real  vertices. 

In  MACH2,  boundary  conditions  throughout  the  code  are  applied  using 
ghost  cells.  The  principal  advantage  of  this  choice  is  that  the  expres¬ 
sions  for  the  exterior  data  are  simpler  to  code  than  the  boundary  differ¬ 
ence  equations.  The  principal  disadvantage  is  the  extra  storage  required. 

A  lesser  disadvantage,  the  extra  computations  required  at  the  boundary,  is 
of  negligible  effect,  since  these  are  carried  out  in  a  vectorizable  loop. 
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Since  <;impler  code  is  rnore  reliable,  the  choice  has  the  effect  of  using 
slightly  greater  computer  resources  to  reduce  the  physicist  resources 
required  to  run  a  given  simulation. 

Ghost  cell  boundary  condition  application  works  well  with  finite 
volume  differencing,  because  only  first  order  differences  are  computed.  To 
correctly  compute  vertex  differences  of  cell  centered  quantities  along  the 
boundary,  it  is  necessary  to  fill  that  boundary's  ghost  cells  with  appro¬ 
priate  data  before  the  difference  computation.  However,  the  cell  centered 
differences  of  vertex  data  along  a  boundary  do  not  reference  ghost  vertex 
data;  thus  it  is  not  usually  necessary  to  fill  ghost  vertices  with  data 
before  performing  the  differences. 

9.1  LOCATION  OF  GHOST  VERTICES 


The  geometric  coefficients  described  in  Section  7.4  are  the  principal 
exception  to  the  rule  just  expressed.  They  are  differences  of  the  vertex 
coordinates,  and  are  thus  cell  quantities.  Their  values  are  required  in 
the  ghost  cells  to  compute  boundary  differences  such  as  the  gradient  or 
curl  of  cell  centered  quantities,  as  in  Equations  121  and  133.  While  these 
ghost  cell  values  could  be  determined  by  symmetry  considerations,  it  is 
more  reliable  to  set  ghost  vertex  locations  and  compute  the  ghost  cell 
values  of  the  geometric  quantities  exactly  as  for  the  real  cells. 


9.1.1  Boundary  Ghost  Vertex  Location--There  are  two  cases  for  the 
formation  of  ghost  vertices  at  a  given  boundary  of  a  given  block.  In  the 
first,  another  block  abuts  this  one  across  this  boundary.  The  relevant 
physical  boundary  condition  is  continuity  of  physical  data,  so  the  ghost 
vertex  positions  in  this  block  are  copied  from  the  first  interior  row  of 
vertices  in  the  adjoining  block.  Thus,  the  ghost  cells  exactly  overlay  the 


first  row  of  cells  in  the  neighbor  block  as  in  Figure  15.  The  proper  cou¬ 
pling  of  physical  processes  is  ensured  by  the  simple  expedient  of  copying 
data  from  the  neighbor  block  into  the  ghost  cells  of  the  present  block. 


In  the  second  case,  i.e.,  no  adjoining  block  across  this  boundary,  the 
ghost  vertex  locations  are  created  using  the  locations  of  the  vertices  in 
this  block.  Since  most  of  the  physical  boundary  conditions  described  in 
Section  3  involve  normal  and  tangential  derivatives,  the  ghost  vertices  are 
positioned  to  facilitate  the  computation  of  those.  The  position  of  each 
vertex  in  the  first  interior  row  is  reflected  out  through  the  tangent  to 
the  boundary  at  the  nearest  boundary  vertex  to  form  the  ghost  vertex  at 
that  row  or  column  (see  Figure  16).  The  tangent  vector  used  is  the 
centered  difference  of  the  adjacent  boundary  vertices  after  being  rotated 
outward  and  normalized  to  unit  length.  At  the  ends  of  the  boundary,  the 
one  sided  difference  is  used  for  the  tangent  instead  of  the  centered 
di fference. 

If  the  boundary  turns  too  sharply,  then  this  process  will  generate 
ghost  cells  which  are  bow  tied  as  in  Figure  17.  In  most  instances,  this 
will  be  avoided  if  the  distance  from  the  boundary  to  the  next  gridline 
inward  is  less  than  the  radius  of  curvature  of  the  boundary. 

9.1.2  Corner  Ghost  Vertex  Location--The  scheme  described  above  does 
not  determine  the  ghost  corner  vertex  location.  It  is  important  to 
position  it  so  that  the  corner  ghost  cell  formed  should,  when  possible, 
overlay  physically  corresponding  cells  in  neighbor  blocks,  be  they  ghost  or 
real.  There  are  four  cases,  corresponding  to  the  number  of  blocks  which 
may  come  together  at  a  corner. 

If  there  is  only  one  block  incident  at  the  given  corner,  the  corner 
ghost  vertex  is  positioned  so  as  to  make  the  corner  ghost  cell  a  parallelo¬ 
gram  as  shown  in  Figure  18a.  If  the  interior  angle  of  the  block  corner  is 
less  than  60  degrees,  such  as  in  Figure  18b,  the  resulting  corner  ghost 
cell  will  have  negative  area  and  volume.  This  restricts  the  range  of 
probl ems . 

When  two  blocks  meet  at  a  corner,  they  are  required  to  also  meet  in  a 
whole  edge.  Then  the  ghost  corner  vertex  position  is  chosen  to  be  the 
same  as  the  position  of  the  ghost  vertex  from  the  neighbor  block  adjacent 
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16.  Ghost  cells  on  boundary  with  no  neighbor  block 
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Figure  18.  Corner  ghost  cells  at  one  block  corners  with  interior  angle 
a)  less  than  60  degrees,  b)  greater  than  60  degrees. 
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Figure  20.  Corner  ghost  cells  at  three  block  corner. 
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Figure  22.  Corner  ghost  cells  at  two  block  corner  in  curved  boundary, 


to  the  real  vertex  which  is  itself  adjacent  to  the  real  corner  vertex  along 
the  non-coincident  edge  as  in  Figure  19.  Note  that  this  means  that  all 
boundary  ghost  vertices  in  all  blocks  must  be  positioned  before  any  corner 
ghost  vertex  positions  are  determined. 

When  three  blocks  come  together  at  a  corner,  the  effort  required  to 
make  the  ghost  vertices  overlap  near  that  corner  is  far  too  great,  even  to 
achieve  validity  only  for  a  limited  set  of  grids.  Hence,  the  ghost  vertex 
in  each  block  is  positioned,  just  as  for  the  case  of  a  one  block  corner,  by 
making  the  ghost  corner  cell  a  parallelogram.  The  four  cells  adjacent  to 
the  corner  may  thus  be  quite  different  from  block  to  block  (see  Figure  20). 
This  results  in  different  values  for  the  geometric  coefficients  in  cells 
near  the  corner,  and  different  values  of  vertex  differences  such  as  current 
density  at  the  real  vertex.  For  some  quantities,  these  different  values 
may  be  sensible,  as  they  may  be  considered  to  be  the  different  limits  of 
that  physical  quantity  as  the  corner  is  approached  from  different 
directions.  For  others,  such  behavior  may  be  unphysical  or  may  generate 
numerical  instability.  In  those  cases,  it  may  be  appropriate  to  average 
the  different  representations  together  and  set  all  of  them  to  the  same 
value. 

The  case  of  four  blocks  meeting  at  a  corner  is  easily  handled.  The 
ghost  corner  vertex  is  placed  at  the  natural  position,  coincident  with  the 
appropriate  real  vertex  in  the  block  diagonally  opposite  the  one  under  con¬ 
sideration  as  shown  in  Figure  21. 

When  two  blocks  come  together  at  a  corner  and  the  non-coincident  edge 
is  not  a  straight  line,  a  further  problem  exists  with  the  ghost  vertex 
adjacent  to  the  real  corner  vertex  off  the  non-coincident  edge.  Recall 
that  the  tangent  vector  at  the  end  of  a  boundary  points  directly  toward  the 
next  vertex  along  that  boundary  in  that  block.  Hence,  the  tangent  vectors 
at  the  corner  points  of  the  two  blocks  are  not  identical.  As  Figure  22 
shows,  the  two  representations  of  that  ghost  vertex  will  then  have  differ¬ 
ent  positions  due  the  use  of  different  tangent  vectors  to  reflect  the 
interior  vertex  outward.  In  this  case  it  is  advisable  to  average  these  two 
positions  and  use  the  result  for  both. 


Of  course,  there  is  an  effective  radius  of  curvature  at  such  a  corner, 
and  the  ghost  cells  may  bow  tie  there  unless  care  is  taken.  That  radius  is 
approximately 

d  1  +  d  1  f  1  1  1 
2  sTn  0 

where  9  is  the  angular  change  of  the  boundary  at  the  corner,  and  dl  and  dl' 
are  the  distances  from  the  corner  to  the  nearest  neighbor  gridpoints  along 
the  boundary. 

9.2  SCALAR  BOUNDARY  CONDITIONS 

The  physical  boundary  conditions  described  in  Section  3  apply  to  both 
scalar  and  vector  functions.  Most  of  those  applicable  to  scalar  functions 
are  of  one  of  two  forms:  either 


n  •  =  0 


on  3R 


(172) 


on  3R 


(173) 


for  some  unknown  scalar  function  (})  and  some  known  scalar  function  f.  In 
almost  all  cases,  <{i  is  a  cell  centered  quantity. 


Equation  172  may  then  be  implemented  by  simply  setting  the  ghost  cell 
values  of  equal  to  its  values  in  the  adjacent  real  boundary  cells. 

Because  of  the  way  the  ghost  cells  were  constructed  above,  the  normal  deri¬ 
vative  of  at  the  boundary  will  be  zero  to  at  least  second  order  in  dl. 

One  special  case  of  particular  interest  is  Equation  25,  the  conducting  wall 
boundary  condition  for  Bg.  This  is  implemented  by  forming  the  cell 
centered  radii  in  the  ghost  cell  and  the  real  cell,  dividing  former  by  the 
latter,  and  setting  the  ghost  cell  value  of  Bg  to  the  product  of  that  ratio 
and  the  real  cell  value  of  B„. 


»  ) 


m 


In  nost  cases  involving  Equation  173,  the  function  f  is  constant.  For 
cell  centered  it  might  be  argued  that  the  correct  way  to  implement  this 
is  to  set  the  ghost  cell  value  of  ^  so  that  the  average  of  the  ghost  and 
real  cell  values  is  f,  i.e.. 


=  2f  - 


where  the  subscripts  r  and  g  stand  for  real  and  ghost,  respecti vely.  How¬ 
ever,  as  the  result  of  this  boundary  condition  will  surely  be  that 


~  (J) 


it  may  be  implemented  nearly  as  accurately  by  simply  applying 


,g  =  f 


(175) 


(176) 


Once  again  Bg  is  a  special  case.  Equation  26,  one  form  of  Ampere's  Law,  is 
implemented  by  setting 


where  r^  is  the  cell  centered  radius  and  I  is  the  total  current  flowing 
between  this  boundary  and  the  axis  of  cylindrical  symmetry. 


9.3  VECTOR  BOUNDARY  CONDITIONS 

The  vector  boundary  conditions  of  Section  3  are  mostly  of  one  of  the 
following  three  forms 


n  •  V  =  0  on  3R 


(178) 


or 


V  =  0  on  3R 


j  1  •  V  ilj2  •  V )  =  0 


on  dR 


(180) 


where  each  of  ui  and  U2  is  either  n  or  t.  Two  such  conditions  are  applied 
to  any  vector  field  on  each  boundary,  with  Equation  179  counting  as  two. 

9.3.1  Vertex  Centered  Fields--If  V  is  vertex  centered,  then  Equations 
178  and  179  may  be  applied  directly  to  the  real  boundary  vertex  data. 
Equation  178  is  applied  by  replacing  the  values  of  V  using 


V  (t  •  V)  t 


(181) 


No  ghost  vertex  values  need  to  be  set.  The  tangent  T  is  computed  using 
"upwind"  differencing  for  quantities  "V  which  have  direct  effect  on  the  grid 
positioning,  such  as  the  grid  displacement  or  fluid  velocity.  This  means 
that  the  tangent  is  computed  using  a  one  sided  difference  in  the  direction 
of  V,  and  amounts  to  having  each  gridpoint  follow  the  gridpoint  ahead. 

This  is  necessary  because  centered  differencing  in  these  cases  produces  an 
instability  which  destroys  the  smoothness  of  the  grid  on  a  moving  boundary. 
For  quantities  which  are  not  closely  coupled  to  the  grid  position,  centered 
differencing  for  the  tangent  suffices. 

9.3.2  Cell  Centered  Fields — For  cell  centered  quantities  V,  Equations 
178  and  180  are  the  most  commonly  used.  One  condition  on  the  normal  com¬ 
ponent  and  another  on  the  tangential  component  are  required  at  each 
boundary.  This  pair  of  conditions  is  combined  into  a  single  simple  geo¬ 
metric  relationship  between  the  ghost  cell  value  and  the  real  edge  cell 
value  of  the  vector  field.  For  example,  the  perfectly  conducting  wall  con¬ 
dition  for  8^  and  8^,  from  Equation  24,  is 


n  •  B  =  0 
and 


on  3R 


(182) 


n  •  V  (t  •  B)  =  0 

These  nay  be  accomplished  by  setting  the  normal  component  of  (B^,  8^)  in 
the  ghost  cell  equal  to  the  negative  of  its  value  in  the  adjacent  real  edge 
cell,  and  the  tangential  component  equal  to  its  real  edge  value. 

Geometrically,  this  is  equivalent  to  reflecting  (B^,  B^)  out  through 
the  tangent  to  boundary.  It  is  by  this  geometric  description,  tangent 
reflection,  that  the  operation  is  identified  in  the  code.  Other  similar 
geometric  operations  which  are  used  include  normal  reflection,  and  normal 
and  tangent  projection. 

9.4  SEQUENCING  BOUNDARY  CONDITIONS  TO  CONTROL  CORNER  VALUES 

Which  of  the  two  boundary  conditions  to  apply  at  a  corner  is  too 
difficult  a  question  to  answer  a  priori.  Therefore,  this  choice  is  left 
open  so  that  it  may  be  settled  differently  in  different  cases  by  careful 
thought  or  experiment.  The  choice  is  made  for  each  physical  boundary  con¬ 
dition  by  specifying  an  order  for  the  boundaries  of  each  block.  Since,  the 
boundary  conditions  are  set  from  one  ghost  corner  to  another,  the  last  con¬ 
dition  applied  on  the  two  boundaries  meeting  at  a  given  corner  determines 
that  corner's  ghost  cell  value. 
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